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A Wing-Body Problem in a Supersonic 
Conical Flow 


S. H. BROWNE,* L. FRIEDMAN,* anv I. HODES* 


North American Aviation, Fne. 


SUMMARY 


The linearized potential equation for steady-state supersonic 
tiow is solved for a slender conical body with a thin symmetrical 
wedge wing extending to infinity in the direction of flow. The 
potential function is divided into symmetric and antisymmetric 
parts, the first of which corresponds to a flow arising from thick- 
ness distribution with no lift, while the second arises from a super- 
posed cross-flow that determines the angle of attack and gives 
rise to lift. Pressure distributions are computed for both cases. 
In the symmetric case the resultant pressure is less than the sum of 
wing and body contributions but never less than that of either 
element alone. In the antisymmetric case, if the wing is wide, 
the pressure and lift curve slope are nearly those of a ‘wing of 
the same plan form and angle of attack. As the wing is 
made narrower, the effective angle of attack of the system 
increases. 


(1) INTRODUCTION 


bee PARTIAL DIFFERENTIAL EQUATION satisfied by 
the velocity potential for steady-state compres- 
sible irrotational flow may be made linear by the as- 
sumption that there are everywhere only small per- 
turbations in the fundamental constant velocity flow. 
Consequently, in this theory it is possible to build up 
the solution for flow around a given body or system of 
bodies by superposition, of solutions corresponding to 
simpler systems. However, if the linearized solution 
for a body of revolution is added to that for a wing, the 
velocity field of each will disturb that of the other, and 
the velocity field of the combination will not satisfy 
the boundary conditions for the physical system of the 
body with a wing. To obviate this difficulty, other 
solutions of the differential equation are added such 
that the total potential and its derivatives assume the 
required values everywhere on the boundaries. These 
additional solutions will be collectively termed the 
interaction potential. 


Presented at the Aerodynamics Session, Sixteenth Annual 
Meeting, I.A.S., New York, January 26-29, 1948, 
* Aerophysics Laboratory. 


In general, the problem of determining the inter- 
action for any given combination of bodies is somewhat 
complicated. As a first step, however, much useful 
information may be obtained from a treatment of the 
special case of a system for which the resultant flow is 
conical. To this end, the system selected for study was 
a slender conical body with a thin symmetrical wedge 
wing, both of which extend to infinity in the direction 
of flow (Fig. 1).* The velocity field for each of these 
alone in a supersonic flow is known from the work of a 
number of investigators.'~* It remains, therefore, to 
determine the interaction of the combination. 


(2) GENERAL SOLUTION FOR CONICAL FLow 


In cylindrical coordinates, the linearized potential 
equation for supersonic flow is 


1 1 
+ br + — (M* — 1) bn = 0 (1) 


where ¢ is the potential of the perturbation velocity and 


L MACH CONE 


Fic. 1. Conical wing-body system. 
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MACH CONE 
Fic. 2. Coordinate system. 


M is the Mach Number of the undisturbed flow, the 
direction of which is taken parallel to the z-axis, *In a 
recent paper by Hayes‘ it is shown that with the intro- 
duction of a conical (nonorthogonal) coordinate system, 
(Fig. 2) through the relation, 


t =rVM? — 1/z (2) 


solutions of Eq. (1) may be put iato a form that is 
especially adapted to the description of conical flow 
fields. Under this transformation, Eq. (1) becomes 


1 1 
(1 — #)bu + (1 — 2#)¢, + + 
— = (3) 
By separation of variables, solutions of the form 
are found, where 7) satisfies the equation 


[1 + 2(m — — 
n(n — 1)?)T=0 (5) 


Solutions of Eq. (5) for integral values of m and non- 
negative integral values of m are discussed by Hayes, 


who found that the functions T,,,,(¢) may be expressed . 


in terms of hypergeometric functions. They fall into 
two classes, designated as Class I and Class II, depend- 
ing upon their behavior at ¢ = 1 (i.e., on the Mach cone 
arising from a disturbance at the origin). Class I 


functions, which will be denoted by S,,,(¢) in this 
paper, exist everywhere in space, while Class II func- 
tions, here designated R,,(¢), are real only for t < 1 
and vanish at ¢ = 1, corresponding to the physical be- 
havior of small disturbances in a supersonic flow. For 
m > |n + 1/,|, the functions may be expressed in 
terms of terminating polynomials in /*; 


for m < 


1948 


\n + 1/,|, they may be logarithmic and are directly re- 
lated to the associated Legendre functions, The gener- 
ating equations deduced by Hayes may be modified to 
give recurrence relations by utilizing certain well- 
known properties of the hypergeometric functions. 
These are { 


(2n + 1)Ta-1, m = "Tn, m — U(dTy, m/dt) (6a) 
(2n + 1)Tp-1,m41 = (m/t)Tam — (dT am/dt) (6b) 
(2n + —(m/t)Tam (dT (6c) 


The perturbation velocity is the gradient of the per- 
turbation potential. If u, v, and w are the x-, y-, and 
z-components of this velocity and g, and g, are the 
radial and tangential components, respectively, then, 
in terms of the conical coordinates used, 


qr = (V M? — 1/2)(0¢/Ot) (7a) 

go = (V M? — 1/tz)(0g/08) (7b) 
and 

w = (0¢/dz) — (t/z)(0¢/dt) (7c) 


In a conical flow, velocity components are constant 
along rays passing through the origin. As can be seen 
from Eqs. (4) and (7a) to (7c), the potential for a conical 
field may be expressed by appropriate superposition of 
solutions of Eq. (3) for which m = 1. Since any deriva- 
tive of ¢ with respect to a rectangular Cartesian co- 
ordinate also satisfies Eq. (1) or (3), Cartesian com- 
ponents of velocity in a conical field may be described 
by similar superpositions with » = 0. In these cases 
especially simple relations exist between the R,zm and 
Snm functions. Some of these relations, together with 
a table of the first few Ri, and Ro, functions are given 
in the Appendix. 

In the case of conical flows, an alternate method of 
formulation, discovered by Busemann,® is_ useful. 
Under the change of variables 


t Qs 


Eq. (3), with m = 0, reduces to Laplace’s equation in the 
polar coordinates sand 6. Thus, for conical flow, in the 
region inside the Mach cone (t < 1), the velocity com- 
ponents u, v, and w individually satisfy Laplace’s equa- 
tion in these variables and may be represented as the 
real or imaginary part of an analytic function of the 
complex variable « = se’. In particular, if 


(8) 


or 


w = (1/VM? — 1)Re fie) (9a) 
it can be shown®’ that 

u = —Ref [(1 + ¢)/2]df (9b) 
and 

v= —Reif[(1 — (9c) 


with a suitably chosen constant of integration. . 
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and the Mach cone as the unit circle. 


A WING-BODY PROBLEM 


The pressure difference from the pressure in the 
undisturbed flow, expressed in terms of the fluid veloci- 
ties, is 

‘p= —pl Vw + [(u? + v*)/2]} 


where p and V are the free-stream density and velocity, 
respectively. When the squares of the components of 
the cross-velocity, u and v, are small compared to Vw, as 
is the case for planar systems, this reduces to 


p= 


For nonplanar bodies, which lie everywhere near the 
z-axis, it may be necessary to use Eq. (10a) to compute 
the pressure. 


(10a) 


—pVw (10b) 


(3) Bounpary ConDITIONS 


The selection of suitable solutions of the type of Eq. 
(4) depends upon the boundary conditions to be satis- 
fied. In general, the condition at the surface of a body 
in a flow is that the resultant fluid velocity vector 
should lie in the tangent plane at each point—that is, 
qn = 0, where g is the total velocity and 7 is a unit 
vector normal to the surface. 


In a conical flow, velocity and pressure are functions 
of tand @ only. A plane in which ¢ and 6 are taken as 
polar variables may be regarded as any downstream 
section of the system normal to the z-axis. In this 
plane, the conical body appears as a circle of radius f 
The leading 
edges of the wing are represented by points on the 
x-axis at +t, and since the wing constitutes a planar 
system, the boundary conditions that exist at its sur- 
faces will be satisfied on the segments of the x-axis 
corresponding to ft < lt < 4. This is equivalent to 
splitting the cone on a diameter and separating the 
two halves to pass the wing through, in contradistinc- 
tion to’ the alternate procedure of removing a portion of 
the cone to permit insertion of the wing (Fig. 3). The 
distinction is not of importance if the semiwedge angle, 
8, is small compared to the semicone angle, ¥, but does 
permit the application of the results to a system con- 
sisting of a split cone and a wing even where 8 is of the 
order of Yo, provided the wing is wide enough. If the 
wing is narrow (4, approaching f) with 6 nearly as 
large as Yo, no part of the system may be considered as 
planar, and a somewhat different treatment would be 
required. 

It is convenient for analysis to separate the potential 
function into symmetric and antisymmetric parts with 
respect to the xz-plane. When the boundary conditions 


at the surfaces of the system also may be expressed as 
a sum of symmetric and antisymmetric terms, this 
procedure separates the general problem into symmetric 
and antisymmetric problems, each of which may be 
solved independently. It is easily shown that the ve- 
locity components u and w have the same type of sym- 
metry as ¢, while that of v is of the opposite type. The 
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symmetric part of ¢ corresponds to a flow arising from 
thickness distribution with no lift distribution,’ while 
the antisymmetric terms represent a superposed cross- 
flow that determines the angle of attack and gives rise 
to lift. 

While the velocity fields for the two types of flow may 
be determined independently, this is not true, in gen- 
eral, of the pressure. For planar systems, where Eq. 
(10b) is used, the pressure has the same symmetry 
properties as w, and the fields due to thickness distribu- 
tion and lift distribution are independent and non- 
interacting. For nonplanar bodies, such as the conical 
body in the system considered here, Eq. (10a) must 
be used to compute the pressure, and interaction be- 
tween the symmetric and antisymmetric velocity fields 
may arise from the quadratic terms in this equation. 

The symmetric problem is simpler and will be dis- 
cussed first. To the approximation that the slopes of 
all surfaces of the system with respect to the direction 
of flow are small, the general boundary condition stated 
above reduces to g, = Vy on the body andv = +V8 
on the wing. Since v is antisymmetric with respect 
to y, it must vanish everywhere on the x-axis outside of 
the system. The boundary condition at the Mach cone 
depends upon the position of the leading edge. If this 
edge lies inside the Mach cone, all velocity components 
vanish for ¢ = 1. If the wing penetrates the Mach 
cone, the velocities must reduce to constant values 
along definite segments of the arc of the unit circle. 
In the latter case, the boundary conditions on the unit 
circle are satisfied by the wing solutions, which will be 
taken as necessary from the literature, and hence need 
not be discussed here. The interaction will always be 
required to vanish for ¢ = 1. However, because of 
differences in the analytical expressions for the velocity 
field of the wing in these two cases, they are treated 
separately. 

In the antisymmetric problem, the wing is considered 
to be a flat plate of zero thickness, coincident with the 
xz-plane. The entire system is placed at an angle of 
attack, a, by superposing a uniform upward component 
of velocity of magnitude V sin a, or Va in linearized 
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theory. In order to give a resultant free-stream velocity 
of V and to make p = 0 in the undisturbed flow when 
Eq. (10a) is used, it is necessary to superpose also a 
uniform axial velocity — V(1 — cos a), approximated 
as —!/2Va?. However, since this term is of second 
order, it need not be considered when computing the 
pressure by Eq. (10b) for planar systems. The bound- 
ary conditions for the antisymmetric flow are, then, 
gr = O for t = ty (since the condition g, = Vy is already 
satisfied by the symmetric part of the total g,), v = 0 
on the x-axis for fh < {t| < t, and v = Va (or q, = 
Va sin 6) for t = 1. When the leading edges of the 
wing lie outside of the Mach cone there are no disturb- 
ances ahead of the wing, so that all velocity compo- 
nents vanish on the x-axis for It| >t, When the system 
is entirely within the Mach cone, however, the value of 
v is not knowh on the x-axis between the Mach cone 
and the leading edges of the wing. Here, since the 
pressure is antisymmetric with respect to y, a condition 
is imposed on w which must vanish on the x-axis for 
th<|t|<1. 

From the above considerations of symmetry and 
leading-edge position, the general problem falls natu- 
rally into two symmetric and two antisymmetric cases, 
each of which will be discussed in the following para- 


graphs. 
(4) SyMMETRIC PROBLEM 


Both cases of the symmetric problem may be solved 
formally as follows. Let the total potential be ex- 


pressed as 
= + dwt or (11) 


where ¢, is the potential for a cone alone, ¢w is that 
for the wing alone, and ¢, is the interaction. The first 
two terms are known, and conditions are sought which 
will permit determination of ¢;. From Eqs. (7a) and 
(7b), the radial and tangential velocities are 
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For 
the: 


= Vie = 1 Odw 


and 
and 


de 


VM? — 1 (2 Odw 2) 


For 6 = 0, g = v = V8 (on the upper surface of the 
wing) for t < |t| < and g = 0 < But this 
conditionis satisfied by the term M? — 1/tz) (O¢w/d8)p PY 
Furthermore, from symmetry, ¢, is not a function of oan 
6. Hence, ¢; must satisfy the condition 


from 


= 0 (18a) 
On the cone, at ¢ = f, the condition g, = Vy is satis. By di 

fied by the first term in Eq. (12b). Conseqtiently, the 

second required condition is 

= (13b) 
From previous considerations, ¢,; may be expressed asa acai 
superposition of R type solutions, - Whe 

=2 (Ansin m0 + B, cos (t) (14 
where 

Substitution of Eq. (14) into Eqs. (13a) and (13b) 

yields A,, = 0 for all m, and 

Cn cos = — Vita i (15) 

where 
dR, and 

Cm Bm ( dt 
and lih<< 
expressic 
_ VM? = 
= ot 


For a given q,,,, the C,, and, hence, the By, may be de- 
termined by the standard technique of Fourier analysis 'rom wh: 


for comy 
wally, | 
Case I. Leading Edges of Wing Outside of Mach Cone 
To calculate g,,, for each wing configuration it is most convenient to use Eqs. (9a), (9b), and (9c), since the 
quired f(e) is already known.’ An f(e) that satisfies the boundary conditions for a wing that penetrates the M 
cone is Case TIT. 
Ve th + 1)? When t 
fwle) = lg 2 Mach cor 
V ty? —1 (e )(e + 1) determine 
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For this case, since the assumptions of the linearized theory require that t << 1, it is sufficiently accurate to use 


the approximate expressions 


Uw = (2V8/z)(t/ti) cos 0 


and 


tw = VB — (2VB/r)(t/t) sin 6 


from which 


Grw = VB sin + (2VB/x)(t/ti) cos 20 
By substituting this value of g-w in Eq. (15), the constants B,, may be found, and the potential ¢, is determined 


from Eq. (14) to be 
4VB 


to 


1 


Ro — 3) cos 20 + 75 Rul?) cos40 +... 


Wr= 
M? — 1 


Case II. Leading Edges of Wing Inside Mach Cone 


4V, t 


VB 


6 


where 5; is the value of s in the e plane corresponding to ¢ = ¢; in the ¢, 6 plane. 


1 — t2V1 — cos 6 


by differentiation and application of Eq. (6), it is found that the axial perturbation velocity due to the interaction is 


3 (* 
2° 


1 0° 
Rez (t) cos 26 + cos 40 +.. ] 


- When the wing lies entirely within the Mach cone, a suitable f(e) is 
= (VBh/eV1 — log — 51°)/(1 — 


From this equation, 


uy = = 


and g,,, may be found from the relation 


for computations. 
ically. In this way the curves of Fig. 7 are plotted. 


t? + t,? + cos? 6) 


, sin ) 


t,? t?t,? sin? 6 


= Uw cos + vw sin 


lif << 4 (leading edge of wing near Mach cone), these velocity components are given by the same approximate 
expressions as in Case I, and the same equation results for wy. 


Ce (2/r) fo" — [(Qryp) t= M? 1] cos m6 dé 
Co = (1/4) So" — M? — 1) 


from which B,, and By follow. These integrals may be evaluated exactly, but the resulting expressions are not useful 
In any particular application it is more convenient to perform the Fourier analysis numer- 


More generally, however, 


(5) ANTISYMMETRIC PROBLEM 


Case III. Leading Edges of Wing Outside of Mach Cone 


When the leading edges of the wing lie outside of the 
Mach cone, the pressure due to angle of attack may be 
determined by the method previously used. As was 
toted above, the potential function for this problem 
ust be antisymmetric with respect to the xz-plane. 
However, the antisymmetric problem may be replaced 
bya symmetric problem that can be solved with less 
ificulty. From the assumption that the flow is irro- 
lational, the condition v = constant may be replaced by 
the equivalent condition 0Ow/dn = 0, which is a natural 
ymmetry condition on w at a free boundary. Thus, 


if the wing is replaced by a free boundary and a solu- 
tion found which is symmetric in y, it is the correct 
solution for positive y, and the proper solution for nega- 
tive y is obtained by appropriate changes in sign. 

In formulating the potential functions for this prob- 
lem, it is more convenient not to include explicitly in 
the total potential the term corresponding to the cone 
alone, since the boundary condition on v, determined 
by angle of attack, is satisfied by the wing solution. 
Thus, 


(16) 


where ¢* has the same form as ¢, in Eq. (14). It will 
be’convenient, however, to refer to ¢* as an interaction 


> = ow + ¢* 


S PLANE 
Fic. 4. 


potential. Application of the boundary conditions 
for this case yields Eqs. (13a) and (13b) with @; re- 
placed by ¢*. From the first of these, again, A, = 0 
for all m, while from the second, Eq. (15) is obtained, 
as before, and the constants B,, will be determined by a 
Fourier series expansion of the appropriate velocity 
function. Since the functions R,,,, and their derivatives 
vanish identically at ¢ = 1, there is no contribution 
from ¢* to the upwash at the Mach cone. Thus the 
angle of attack of the total system is the same as that 
corresponding to the wing alone. 

To calculate g,,, for this case, the fw(e) of Case I, 
with 6 replaced by —a, may be used. Eqs. (9b) and 
(9c) now lead to 


sin 6 cos 6 


t 

Vi-# 
2 


sin (6 — @) sin (6 + 8) 


and 


t 


1—# 


Va 
Uw = — tan-! 
T 


sin (6 — 6) sin (6 + 6) 


The difference in the form of vw in this case as com- 
pared with that of Case I arises from the difference in 
boundary conditions on vy in the two cases. Since t, > 
1, the approximate expressions 

Uw = cos 
and 
tw = (2Va/r) (t/t) sin 6 
are sufficiently accurate fort<< 1. Then 


Qry = —(2Va/r) (t/t) cos 20 


and the resulting interaction potential is 
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4 
h 


VM? —1 


From this ¢* 


3Va to! 
t) cos 26 
rVM?—-14 


Thus, in this case, the interaction is negligibly small, 
and the velocity field of the system is essentially that 
due to the wing alone. 


Case IV. Leading Edges of Wing Inside of Mach Cone 


When the system lies entirely within the Mach cone, 
the boundary condition w = 0 on the x-axis for th; < |t|< 
1 makes difficult the selection of a suitable velocity 
potential. The use of the formulation given by 
Eqs. (9) leads more directly to a solution in this case. 
Let the f(e) for the system be expressed as 


= fale) + f*() (17) 


where fw(e) is the function for the wing at an angle of 
attack y. The value of y will not yet be specified, but 
it should be noted that y may differ from a, the geo- 
metric angle of attack of the system. 

The function f(e) is conveniently expressed in the 


¢-plane where 
1 


Note that on the real axis § = ¢. Fig. 4 shows the 
system in the ¢-plane. The body trace remains circular 
of radius t = f) and with o = @ under the assumption 
th << 1. The trace of the Mach cone includes all of 
the real axis except the segment through the origin 
extending from —1 to +1. 

The term fyw(e) in Eq. (17) is known,’ 


fale) = 
— 42) V1 (7/42) 


where E(V1 — t) is the complete elliptic integral 
of the second kind. It is also known that this function 
leads to tw = Vy on the Mach cone, where vp has been 
set equal to zero on the wing. It is now possible to 
determine the boundary conditions on f*(e). It must 
give rise to the following relations: 


(18) 


(19) 


(a) On the real axis for {¢| 2 1,u* = w* = 0. 

(b) On the real axis for 4, < |¢| < 1, w* = 0. 

(c) On the real axis for < < th, 0w*/On = 0. 
(d) On the real axis for 2 1,v* = Via — ¥). 
(e) On the body, for |t|= = —Gry- 


To construct an f*(e) that satisfies these conditions, 
it is convenient to define the T7,,,, functions in the com- 
plex plane. Since ¢ = ¢ on the real axis, from the 
principle of analytic continuation T,m({) is obtained 
from T,»(t) by replacement of ¢ by ¢. A number of 
useful relations involving the Ro,(¢) and Som(f) may 
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then be derived. Some of these are given in the Ap- 
pendix. The boundary conditions (a) through (c) 
above are satisfied if f*(e) is taken as 


= BomRom(/t:) 


where the Bo, are real constants. This may be veri- 


’ fied by inspection if Eqs. (9b) and (9c) are rewritten in 


the form 
(21a) 


(21b) 


u = df 
v= —Reif (V1 — df 
and the normal Sitesi of VM? — 1 w* on the real 
axis is written as 
(0/0n)Re f*(e) = Re i(d/df)f*(e) 


Conditions (d) and (e) may now be used to determine 
the constants Bo, and y. For this purpose, it is neces- 


+ The series 2Bo,Ro», for | t/t| <1 has been proved both ab- 
solutely and uniformly convergent. See Eq. (24d) below. 


= Via — y)sin@ — {Re Bom 


To apply condition (e), an expression for g,,, near the 


Eqs. (21a) and (21b) from the origin under the assumption |¢{| << 1. 


nomial theorem, 


Vy 1 ( 1:3 ( ] 
+3 cos 4 h cos 64 + 
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sary to evaluate the integrals in Eqs. (21a) and (21b), 
which, for u* and v* near the body, lead to incomplete 
elliptic integrals. The computations may be simpli- 
fied by taking in each case a path of integration from a 
point on the real axis for which ¢ > 1 (point on the 
Mach cone) to one on the real axis near the origin 
(point on the wing), then to a point near the origin off 
the real axis. The first portion of this path then con- 
tributes functions of complete elliptic integrals. To 
evaluate the integrals along the second portion of the 
path, it is helpful to make use of the approximate 
relation 


dRom(§/t1) = ty"dRom(¢) 


correct for m = 2 to within 4 per cent for values of 


t/t, up to 0.5 and within 9 per cent for t/t; = 0.75. 
For this integration, then 
df*(e) = 
and the resuiting q,* is 
dRip(t) ‘ 


origin is required. This may be found by integrating 
The result is, after expansion by the bi- 


(23) 


To the approximation desired, only the first two terms need be retained. Comparison of Eqs. (22) and (23), with 


t= to, yields the relations 


Boz = —(3/2)(to/t)*( (24a) 
Bog = Vaty/AEp) (24b) 
Bos = — (21/256) (to/ti)'*( Vaty/AEo) (24c) 
Bo, = —{(2k + 1)(2k — 2)1/2%-*[(k — 1)!]?} (24d) 
y = a/A (24e) 
where 
and where Jo, J4, and J, are 
1 dt 1— + 15 1 o( 
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Here Ey(V 1 — t?) and Ko( eT. t,") are complete elliptic integrals of the first and second kind, respectively. It : 
should be noted that the values of J2, J, and J, depend only on /; and not on f. 
From Eq. (9a), using these results, the axial component of the total perturbation velocity valid for t/t; < 0.75is; " 
(a) On the wing 
VM? — 1 V1 — 2\n) \e a(t) (ek - 
sali) (8 
16— — sot 
i 
(b) On the body 0s 
AE) VM? — 1 cos 260 + - ( cos 46) + 16 (25h) 
In the case of t/t; — 1, the system is small enough W = Va (4 + h4) — 2t*t,? that ¢ 
compared to the Mach cone, so that the latter may be _ VM? — 1 Lev 2(4,2 — 42)? — 4208 + = | of thi 
considered infinitely far away. This method of treat- comp 
ing low aspect ratio systems, due to Jones,* reduces (28) free-st 
the lift problem to one of an incompressible cross flow havior 
of magnitude Va. In such a flow, the complex velocity (b) on the body, wing i 
potential ¢ is defined by The 
W = Va : wide-v 
do/de = u — w V (h? + — cos? @ those 
the res 
Using Eqs. (9a) and (9b), with «? << 1, yields 
de = —! 
For purposes of discussion, the results of the above] Thus, 
Solving for f and then integrating by parts, analysis have been applied to calculate the pressure dis-] cone), 
fl). = —2(do/de) + 26 (26) tributions for particular examples of each of the four 
. cases treated. For each case, an appropriate 4 anda 
Under the Joukowski transformation to = 0.1322 (corresponding to a cone of 10° semiverter} In C 
£ = Wale + (s02/e)] angle at M = 1.25) are selected. The Pressure dis ing to 
tributions that arise in the two symmetric, wide-wing] metric 
the wing body system transforms into a flat plate, in cases are shown in Figs. 5 and 6, where pressure co} is show 
efficient, Cp = p/(p/2)V’, at zero angle of attack, is tion cal 


the — plane, whose potential is given by 
= UVaV# — a? (27) 


where a (s1? + 50?)/2s,. Eqs. (26), (27), and (9a) 
yield: (a) on the wing, 


plotted as a function of ¢. 


The contribution of the plotted 


interaction is of the same form in each case. It can bef of these 
seen that the resultant pressure is less than the sum off the pre: 


the wing and body contributions but is never less than 
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that of either element alone. Thus the drag of a system 
of this type would be greater than the sum of the drags 
computed for each element considered alone in the given 
free-stream flow. As can be seen from Fig. 7, the be- 
havior of the system is essentially unchanged as the 
wing is made narrower. 

The antisymmetric pressure distributions for both 
wide-wing cases are to within a close approximation 
those that occur on a wing alone (extending through 
the region occupied by the body) at the angle of attack 
of the system. The lift curve slope, dC,/da, to the 
same degree of approximation, is that of the wing alone. 
Thus, for Case III (leading-edges outside the Mach 


cone), 
4/V M? — 1 


In Case IV, two conditions are distinguished accord- 
ing to the aspect ratio of the wing. The antisym- 
metric pressure distribution on a system f)/t; = 0.661 
is shown in Fig. 8a. In Fig. 8b the pressure distribu- 
tion calculated from the incompressible flow theory is 
plotted for a system with &/t; = 0.9. Comparison 
of these two curves shows a smooth transition toward 
the pressure distribution for a cone alone. 

From Eq. (25) dC,/da for to/t; < 0.75 is found. 


dC,/da 


(2) 2 (*) (*) | (29) 
The first factor in dC,/da is precisely the lift curve 
slope for the wing alone. The remaining two factors 
may be regarded as representing a decrease in the effec- 
tive angle of attack of the system with increase in the 
tatio t/t; as t; decreases from 1. A physical interpre- 
tation is that the interaction contributes an additional 
downwash at the Mach cone which is effectively 
equivalent to a decrease in the angle of attack of the 
system as compared with a wing of span 4. 
From Eq. (28) dC,/da for t/t; > 0.75 is found. 
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Fic. 8a. Antisymmetric case. Polynomial solution: &/t; = 
0.661, M = 1.25, a = 1.72, Fo = 10°. 
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Fic. 8b. Antisymmetric case. Incompressible flow solution: 
bo/t, = 0.9, M = 1.25, a = 1.72, Yo = 10°. 
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POLYNOMIAL SOLUTION 7- Zero Degree S Functions 
Order 
0 1 
1 1/t 
2 (2—#?) /t? 
3 {1 + 301 — #))/# 
4 [1 + 6(1 — #) + (1 — #)*]/#4 
5 
6 


Fic. 9. Ratio of dCz/da for wing-body system to equivalent 
wing dCL/da. t& = 0.13822. M = 1.25. 


[1 + 10(1 — #) + 5(1 — #)4)/8 
[1 + 15(1 — #) + 15(1 — #)?+ (1 — &)*]/8 


da VM?—114/% (1+)? 


Zero Degree R Functions 
+ Sot Th 
2 2 the cc 
2 ty? th? + ty? #[1 + 1/:(1 — 
—4 
Fig. 9 shows how the two methods give curves of dC, /da 
which fair smoothly into each other for t/t; ~ 0.75. 5 — + 201 — +3401 
The lift curve slope for both the cone alone and the 6 to] — [1 + 0/3) — 2) + (1 — #3] ing hi 
wing alone is 27t/V M* — 1, where ¢ represents or 
t, when t << 1. Since for the system the limiting 4syptotic Values deals 
values fo/t; = O and &/t; = 1 represent the cases of The asymptotic values of dR,,,/dt as t — 0 are given been i 
wing alone and cone alone, respectively, the ratio form from0 to 6. 
(dC,/da)/(dCy/da)wing must be 1 at these extreme been a 
values. Fig. 9 exhibits this behavior. dRy/dt = 1/t 
dRu/dt = —(1/t2) — tanh-'V1 — 2 
APPENDIX dR»/dt = —2/t5 
dRy/dt = —3/t* 
The real and imaginary parts of the R and S functions dRy/dt = —24/565 _— 
have an especially simple form. dRy/dt = —8 /ts either 
dRy/dt = —96/7t' terse! 
Rom({) = Rom(t) cos — — Som(t) sin zatior 
ws R ; greate 
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Helicopter Control and Stability in 
Hovering Flight’ 


R. H. MILLERt 
Massachusetts Institute of Technology 


SUMMARY 


The stability and control characteristics of several different 
types of helicopters have been evaluated and compared with the 
characteristics of a conventional airplane. It is concluded that 
some form of automatic stabilization is indicated if more desirable 
handling qualities are to be attained, but, if stability is not to be 
achieved at the expense of controllability, means of increasing 
the control effectiveness is evidently necessary. Methods of ac- 
complishing this are discussed, and it is shown that, in some cases, 
response and stability characteristics approaching those of an air- 
plane may be achieved. The equations of motion for the heli- 
copter and rotor blades have been developed in nondimensional 
form in such a way as to include the effects of flapping and feather- 
ing hinge offsets and restraint about these hinges. Both single 
and twin rotor helicopters have been considered and the complete 
solution simplified, without sacrificing engineering accuracy, 
wherever possible. The effects of automatic stabilization have 
been included in the solution in a manner that, it is believed, 
will cover most applications, and an elementary servo system has 
been analyzed by means of the frequency response method. 


(1) Discussion 


T= NECESSITY OF DEVELOPING MEANS for simplifying 
and improving helicopter handling characteristics 
either by design changes modifying the basic parame- 
ters effecting stability and control or by artificial stabili- 
zation, is becoming increasingly apparent with the 
greater development of helicopter utility. An at- 
tempt is made in this paper to examine the stability 


‘characteristics and control response of helicopters with 


and without automatic stabilization and with several 
different rotor configurations including single and dual 
rotors, rotors with various degrees of flapping and 
feathering hinge offsets, and rigid feathering rotors. 
The complete statement of the problem results in a 
solution that is too cumbersome for routine application, 
and, consequently, at some sacrifice in accuracy, a more 
simplified solution has been obtained which, however, is 
believed to be adequate for design purposes. 

In order to clarify the control problem and to provide 
some basis for comparison, the initial response of a con- 
ventional airplane to abrupt control manipulation will 
be considered first: In this and all subsequent analyses 
the transient response to an instantaneous control de- 
flection (unit step function), followed by instantaneous 
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neutralization of the controls 1 sec. later, will be used 
to illustrate the control and stability characteristics. 
The response when the initial control deflection is main- 


‘ tained beyond 1 sec. is indicated by a broken line. Fig. 


1 shows the change in angle of attack relative to the 
flight path, for unit elevator deflection, for a typical 
airplane flying at 100 m.p.h. It will be noted that the 
angle of attack reaches a definite value and stays there 
if the control deflection is maintained, a characteristic 
of static stability. The initial response is rapid, and 
the oscillations are so heavily damped as to be almost 
unnoticeable. When the controls are neutralized the 
angle of attack starts returning almost instantaneously 
to the trim position and shows no tendency to continue 
in the direction of the initial displacement—in other 
words, there is practically no “‘overshoot.” 


3 


nm 


2 3 
TIME - SECONDS 


Fic. 1. Initial response of airplane following abrupt elevator 
deflection, 5.. Change in angle of attack, a, relative to flight 
path. (See Appendix.) Elevator neutralized after 1 sec. 
----Elevator deflection maintained. 


10 


2 3 
TIME - SECONDS 


Fic. 2. Initial response of airplane following abrupt elevator 
deflection, 4,. Change in attitude, 6, relative to earth axes. (See 
Appendix.) 


453 


| 2 

[Pak 

PTT VIC 

7 PU 

7 
a 


454 JOURNAL OF THE AERONAUTICAL SCIENCES—AUGUST, 1948 


10 
4 
Z 
7 
L 
4 
4 
6-0 
2 


TIME - SECONDS 


Fic. 3. Initial response of airplane following abrupt aileron 
deflection, 6,. Change in angle of bank, ¢. (See Appendix.) 


maintained. 


The response in roll to an abrupt aileron deflection 
shown in Fig. 3 is similar to the change in pitching atti- 
tude relative to earth as shown in Fig. 2 and is charac- 
terized by considerably more overshoot than in the 
case of Fig. 1, since the airplane keeps on rolling for a 
short while after the ailerons are neutralized before 
stabilizing at a definite angle of bank. Eventually, the 
resultant sideslip will cause a spirally stable airplane to 
return to level flight, as, in the case of Fig. 2, the 
phugoid oscillation will, at constant power, eventually 
return the airplane to its original attitude relative to 
earth axes. However, the long-period oscillation and 
convergence do not appreciably affect the initial re- 
sponse and have been neglected for the purposes of this 
illustration. Actually, the importance of spiral sta- 
bility is largely a matter of compromise and depends on 
the type of airplane being considered. In general, an 
airplane that must be flown over long periods or by an 
inexperienced pilot, and for which a high degree of 
maneuverability is not an essential requirement, should 
possess both lateral.and longitudinal stability. 
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TIME - SECONDS 


Fic. 4. Initial response of hovering helicopter with articu- 
lated rotor, = 0.033, following control deflection, 6. Change in 
angle of inclination relative to trim position. No automatic 
stabilization. Controls neutralized after 1 sec. ---- Control 


defiection maintained. 


Ailerons neutralized after 1 sec. ----Aileron deflection 


A conventional helicopter without automatic sta- 
bilization devices has little static stability in hovering 
flight, where static stability may be defined as the exist- 
ence of a restoring force proportional to displacement 
and tending to return the helicopter to its original 


vertical attitude relative to earth axes. Consequently, - 


as might be expected by analogy with a neutrally stable 
airplane, satisfactory control is largely a matter of skill 
and anticipation on the part of the pilot. In addition, 
the absence of any appreciable amount of damping in 
pitch or roll results in continuous oscillations arising 
from what static stability does exist. Theory indicates 
these oscillations to be unstable with a period long 
enough to make their detection difficult but short 
enough to seriously influence the control response. 

In Fig. 4 the control response of a typical dual rotor 
helicopter with slightly offset flapping hinges (3.33 per 
cent of the blade radius) isshown. The change in angle 
of pitch of the helicopter, relative to earth axes, per de- 
gree of cyclic pitch variation, maximum negative at 
90° azimuth, has been used in this and all subsequent 
plots. The initial response is slightly less favorable 
than in the case of Fig. 2, as evidenced by the slower 
rate of increase of the angle of pitch, but it is believed 
that the tendency to oscillate with large amplitude 
about the neutral position and the large overshoot after 
neutralization of the controls are the more serious 
phenomena indicated by this plot. While the period is 
long enough to be handled and would probably not be 
noticeable except as an oversensitive delayed response 
to corrective control manipulation, when combined 
with the overshoot (which is a characteristic of poor 
initial response and lack of stability) these control 
characteristics must necessitate a considerable amount 
of skill and familiarization on the part of the pilot and 
would make noncontact flying a complicated operation. 


It should be noted that the lag in response arises not ° 


from any lag in the response of th2 rotor blades to ap- 
plied cyclic pit :h variations, which the theory of section 
II indicates is almost instantaneous, but because of the 
inherently small control force available with a flapping 
rotor combined with the undamped oscillatory response 
to control manipulation. This control force could be 
increased at will by increasing the amount of cyclic 
pitch per control stick deflection, since, because of the 
rapid response of the blades and the human delay 
involved in building up full control stick deflection, it 
is unlikely that the blades will stall even for cyclic 
pitch variations well above the'stalling angle. How- 
ever, such a control magnification has obvious disad- 
vantages, associated with the high Coriolis forces in the 
plane of rotation induced by excessive flapping, and 
would make control in forward flight difficult. : 

Flight experience has shown that, by moving the 
blade center of gravity aft of the feathering axis, the 
control sensitivity can be appreciably increased. This 
is investigated in detail in section III, where it is shown 
that the primary effect of such a shift is a decrease in 
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TABLE 1 
Roots of Stability Equation for Helicopter with a Single Rotor and Various Ratios of Pitching to Rolling Moments of Inertia 


7, =TIy —0.95 = 0.49% +0.17 = 0.512 +0.13 = 0.72% 
Zz = —1.17 + 0.68: +0.17 = 0.53% +0.04 + 0.79% 
Ie —1.70 = 0.78 +0.17 = 0.54 —0.05 = 0.79% 
I; =0 —1.50 +0.17 = 0.544 —0.12 + 0.73% 
Dual rotor (pitch) —-1.0 +0.18 = 0.64 


damping in pitch and, consequently, an apparent in- 
crease in control sensitivity. Conversely, it is shown 
that a forward shift of the blade c.g. or a rearward shift 
of the aerodynamic center relative to the feathering 
axis provides a convenient means of improving the 
helicopter stability characteristics, although at the ex- 
pense of controllability, by increasing the inherent 
damping in pitch of the helicopter and lengthening the 
period of the lightly damped slow oscillation to the 
point where the helicopter approaches the control char- 
acteristics of the airplane in roll (Fig. 3). 

It is difficult to establish whether the response in 
pitch of a helicopter should be compared with the re- 
sponse in attitude (pitch or roll) referred to earth axes 
of the airplane (Figs. 2 and 3) or the response of the 
angle of attack relative to the flight path (Fig. 1). 
Certainly, if comparison is made between Fig. 4 and 
Fig. 1, the helicopter would appear to possess none of 
the requirements necessary for satisfactory handJing 
characteristics. However, it is believed that, if the 
helicopter is to prove satisfactory for long-period opera- 
tion by inexperienced pilots and for noncontact flying, 
it should return positively and rapidly to its initial 
vertical position upon neutralization of the controls, a 
requirement that necessitates an adequate amount of 
static stability and damping and a response similar to 
that of Fig. 1. 

The roots of the complete stability equation indicate, 
for the configuration of Fig. 4, an unstable oscillation 
with a period of 9 sec. doubling its amplitude in 3 sec. 
In the simplified solution this oscillation has a period of 
10 sec. and doubles its amplitude in 3.5 sec. A single 
rotor helicopter would have two similar oscillations in 
combined pitch and roll (see Table 1), except that, if 
the moment of inertia in roll is small compared to that 
in pitch, one of the oscillations will be neutrally stable 
and may even be lightly damped. 

Although in the interests of mathematical simplicity 
the equations of motion developed in section II have 
been derived without direct reference to the underlying 
physical concepts, it may be of interest before investi- 
gating means of improving the control response to ex- 
amine some of the primary factors influencing the sta- 
bility characteristics. Referring to Fig. 5, the inclina- 
tion of the rotor tip path plane, (;, will, in a conven- 
tional helicopter, tend to follow the inclination of the 
helicopter, a, by virtue of the cyclical pitch variation, 
with reference to their initial plane of rotation, imposed 
on the rotor blades by the inclination of the rotor shaft 
axis. Since the rotor thrust in hovering flight is normal 


to the tip path plane, this results in a forward force for 
nose-down pitching which propels the helicopter for- 
ward. Because of the resultant forward speed, and in 
the process of equalizing rolling moments in forward 
flight by flapping, the tip path plane will incline back- 
ward, through a small angle A (see Fig. 5), relative to 
the shaft axis. If the helicopter center of gravity is 
assumed for simplicity to be along the shaft axis, the 
thrust will then no longer pass through the center of 


‘ gravity, and a restoring moment is produced propor- 


tional to the angle of inclination, since the forward 
velocity is also proportional to the angle of inclination. 
Additional restoring moments due to the forward veloc- 
ity arise from the slight rearward inclination of the 
thrust vector relative to the tip path plane as a result 
of coning, and from the profile and induced drag of the 
rotor blades. These restoring moments are all small 
but are the primary source of the oscillation in pitch. 
Damping arises from inertia forces, due to Coriolis ac- 
celerations proportional to the angular velocity in pitch, 
which induce a: rearward inclination of the tip path 
plane proportional to forward rate of pitching. If the 
flapping hinge is offset, the inertia and aerodynamic 
forces at this hinge resulting from the blade flapping 
will also produce restoring and damping moments pro- 
portional to a; and &. 

Several additional forces and moments exist in the 
case of a flapping rotor; however, their actions are all 
similar to those discussed above. 

Evidently, in order to improve the stability charac- 
teristics, the designer has at his disposal several expedi- 
ents, some of which will be discussed below. 

(a) Reduce the cyclic pitch variation imposed on 
the rotor because of helicopter pitching. This may be 
done by means of a universal hub, which would com- 
pletely eliminate this cyclic pitch variation, or in sev- 
eral other ways. Essentially, the result is to increase 
AB by eliminating the tendency of the rotor to follow 
the helicopter, thereby introducing a considerable 


Fic. 5. Relative inclinations of helicopter and tip path plane 
following a disturbance from equilibrium position. 
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Controls Neutralized After One Second 
—-—— Controls Neutralized After One Half Second 
——— Control Deflection Maintained 


Fic. 6. Effect of decoupling rotor and helicopter on initial 
response of helicopter with § = 0.033 following abrupt control 
deflection. Assumed servo gear ratios: k = 1,r = 0,/ = 0. 


amount of ‘‘pendulum stability.” For purposes of this 
investigation it will be assumed that a control displace- 
ment proportional to helicopter pitch and opposite in 
sense to the cyclic pitch variation imposed by the heli- 
copter is obtained with a conventional servomechanism 
actuated by the signal from a displacement gyro. 
Since the rotor now no longer follows the helicopter, 
the upsetting moment that the pilot applies to the 
helicopter by an initial displacement of the tip path 
plane decreases as the helicopter responds, and a reduc- 
tion in control effectiveness may therefore be expected. 
Also, the poorly damped oscillation, although stable, is 
of higher frequency than before. A typical response 
curve using the parameters, including hinge offset, of 
Fig. 4 is shown in Fig. 6 for a helicopter with a gyro 
stabilized rotor and servo gear ratio of 1, assuming no 
servo time lag. Servo gear ratio is here defined as the 
blade pitch change produced by the servo divided by 
the disturbance that activates the servo. [See Eq. 
(26), section III.] If servo time lag exists or is built 
into the system in order to improve the initial response, 
the displacement gyro alone can no longer stabilize the 
oscillations. Fig. 13, which is discussed more fully in 
section III, indicates that a time lag of 0.3 sec. is suf- 
ficient to make the oscillation once more undamped. If 
the gear ratio is kept at the minimum value necessary 
‘to ensure neutral stability (about 0.05 from Fig. 12), 
the initial response is essentially the same as in Fig. 4. 


ie} 2 3 4 5 
TIME - SECONDS 


Fic. 7. Effect of automatic stabilization on initial response of 
helicopter with € = 0.033 following abrupt control deflection. 
Assumed servo gear ratios: k = 0.15,r = 1,1 = 0. 
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(b) A logical next step is the inclusion of lead in the 
servo system either by differentiating the displacement 
signal, by means of a restrained gyro, or by offsetting 
the blade center of gravity relative to the feathering 
axis. Provision of some rate control in this manner 
permits increasing the helicopter damping to any de- 
sired degree, and a control response approaching that 
of Fig. 1 may then be achieved. In order not to reduce 
the control sensitivity excessively, a careful choice of 
the servo system gear ratios is essential. A gear ratio, 
[see Eq. (26)] of about 0.15 on both rate and displace- 
ment signals, (kr = k = 0.15) with servo time lag com- 
pensated for by such additional rate control as may be 
necessary, gives the response shown in Fig. 7. This 
response appears promising, although the reduced con- 
trol effectiveness may reduce the maneuverability of 
the helicopter to an undesirable level. 

In the above discussion it has been assumed that 
the servo system is incorporated in the control system 
in such a manner as to operate in both stick-fixed and 
stick-free conditions—that is, the pilot controls the 
helicopter not by overriding the automatic pilot but 
by changing the linkage between the servomotors and 
the swash plate, which would be analogous to controlling 
an airplane that was stable, stick fixed and free, by 
means of trim tabs. Whether such a procedure is bet- 
ter than, or even as desirable as, allowing the automatic 
pilot to operate only in the stick-free condition would 
have to be established by flight experience. In the 
final analysis, an aircraft with satisfactory handling 
characteristics is one that appeats satisfactory to the 
pilot flying it, and certainly, in the past, airplanes with 
poor stick-fixed stability have proved satisfactory when 
stick-free stability has been obtained by incorporating 
springs or bob weights in the control system or by ad- 
justing the elevator floating angle by changing the type 
of aerodynamic balance. Since pilots fly mostly by 
feel rather than control displacement, a helicopter that 
indicates stick-free stability by resistance to control 
manipulation proportional to displacement and rate of 
displacement is probably satisfactory, although it is 
possible that, in the case of inexperienced pilots, a com- 
bination of both stick-fixed stability and control feel 
may prove desirable. In the present analysis, the con- 
clusions reached may also be applied to a helicopter 
with stick-free stability only, in which case the stabiliz- 
ing devices will not reduce the control sensitivity as in- 
dicated but, by virtue of ensuring stick-free stability, 
may be designed to exert control forces indicating a de- 
gree of stability as shown in the response curves, pro- 
viding always that the exertion of such forces does not 
change the gyro-reference—i.e., that the servo system 
has a ‘“‘memory”’ with respect to either earth axes or 
flight path. It is not the purpose of this paper to in- 
vestigate the details of various mechanical or electro- 
mechanical devices that may be developed to achieve 
such handling characteristics but simply to present a 
theory adequate for determining the behavior of the 
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overall system incorporating such devices, and it is be- 
lieved that this may be done by assuming a conventional 
servo system with time lags and rate and displacement 
responses. 


The response curves may be used as indications of 
the stick-free transient motion following any disturb- 
ance, other than control manipulations, which inclines 
the tip path plane. 


(c) Another method for improving the stability 
characteristics exactly opposite in concept to (a) above 
would be an attempt to eliminate the rearward inclina- 
tion Af of the tip path plane which is the primary 
source of the unstable oscillations. This leads logically 
to the rigid feathering rotor, which, as may be expected, 
shows somewhat better stability characteristics. A 
similar helicopter to that discussed in (a), but with in- 
finitely rigid rotor blades and no flapping hinges, would 
have a control response as indicated by Fig. 8. The 
roots of the stability equation indicate a slightly un- 
stable oscillation, since rotor drag and coning still con- 
tribute restoring moments even though flapping is ab- 
sent, with a period of 11 sec. doubling its amplitude in 
22 sec. Because of the larger control moments theo- 
retically available with this type of rotor, the initial 
response is excellent and the overshoot small. How- 
ever, except in the case of extremely small rotors, it 
would be impossible in practice to design extremely rigid 
rotor blades, since structural considerations alone 
would probably require blades of such flexibility as to 
result in control responses more nearly approaching 
those of Fig. 4 than of Fig. 6. For small rotors this 
remains as a promising type of rotor construction, par- 
ticularly since rotor vibration in forward flight with a 
three-bladed rotor would be at a low level, in the ab- 
sence of any appreciable third harmonic aerodynamic 
excitation that could be eliminated completely on an 
infinitely rigid rotor if rolling moments were canceled 
by dual rotating rotors rather than by feathering the 
blades. In view of the similarity in response of this 
type of rotor and that considered in (d) below, the ef- 
fects of automatic stabilization will not be examined 
directly for this case. The conclusions are the same as 
for (d). 

(d) For larger rotors in which blade flexibility pre- 
vents approaching the case of rigid feathering rotors, 
a further possibility of improving stability and control 
exists in view of the importance of the flapping hinge 
offset on the stability derivatives. By locating the 
flapping hinge well out along the blade radius, the aero- 
dynamic and inertia shears act at a greater arm, and 
both control and stability are improved. Fig. 9 shows 
the response for a helicopter similar to that of Fig. 4 
but with- the flapping hinge eccentricity increased ten 
times to 33.3 per cent of the blade radius. In Fig. 9 the 
tesponse is given for the combination of swash-plate 
inclinations which will result in an initial inclination of 
the tip path plane in the forward direction only (see 
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Fic. 8. 
blades following abrupt control deflection 42. 
stabilization. 


Initial response of helicopter with infinitely rigid fixed 
No automatic 
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Fic. 9. Initial response of helicopter with offset flapping 


hinges, § = 0.33, to abrupt control deflections, 6; and @. No au- 
tomatic stabilization. 
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F Fic. 10. Effect of automatic stabilization on control response 
offhelicopter with offset flapping hinges. Assumed servo gear 
ratios: k = 0.2,r = 0,1] = 0. 


section III). The unstable oscillation has for this case 
a period of 11 sec. doubling its amplitude in 21 sec. 

Of considerable interest is the control response with a 
displacement gyro only operating through a servo sys- 
tem with a gear ratio of 0.2. This response is shown in 
Fig. 10 and closely approaches that of Fig. 1. Time 
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delays are relatively unimportant, Fig. 14 indicating 
that a 3-sec. delay would be necessary before the system 
again becomes unstable. 

Although this type of rotor presents several vibra- 
tion problems associated with the offset flapping hinges, 
both it and the rigid rotor discussed in (a) above appear 
to offer several advantages when a high degree of ma- 
neuverability is required, and it is desired to avoid the 
complication of introducing lead, or rate damping, in 
the automatic servo system. 

(e) Ifitis further desired to eliminate the automatic 
pilot completely, the helicopter stability characteris- 
tics may be improved, and a certain amount of feel in- 
troduced into the control system, by an adjustment of 
the blade center of gravity, or the aerodynamic center, 
relative to the feathering axis. The resultant adverse 
steady-state moments in hovering and forward+flight 
may be balanced out by suitable interblade linkages 
and by careful adjustment of the blade camber, aero- 
dynamic axis, and gravity axis relative to each other, 
as discussed in section ITI. 

It is possible by this method to increase the period of 
the lightly damped oscillation to the point where the 
helicopter control characteristics approach those of the 
airplane in roll (Fig. 3), although the amount of damp- 
ing in pitch necessary to accomplish this will result in a 
considerable reduction in control effectiveness which, 
- however, may in turn be improved by an appreciable 
offset of the flapping hinge. 

The lengthening of the lightly damped oscillation has 
the disadvantage of eliminating the tendency of the 
helicopter to return rapidly to its initial vertical posi- 
tion in hovering—a tendency that may not be essential 
but is probably desirable and which would appear at- 
tainable only by means of a displacement gyro or its 
equivalent, with trim control for forward flight. 


NOMENCLATURE 

A, B,C, D = coefficients defined in Eq. (12) 

E, li’. J, K = elements of coefficients of stability, Eq. (19) 

L, N, P, C = coefficients defined in Eqs. (14) and (15) 

Au to A33 

and 

Ay, to A3; = elements of the characteristic determinant 

C, to Cs = coefficients of stability Eq. (18) 

Coo = profile drag coefficient of blade element 

, ae = moment coefficient of blade element 

H,toH, = coefficients defined in Eqs. (14a) and (15a) 

1s Pp: & = horizontal forces in reference plane acting on 
the rotor, Ibs. 

Ty, to I;3 - = imaginary parts of 4, to Ag; 

Fiods = moments of inertia of helicopter in pitch-and 
roll, Ibs. ft. sec.? 

I = (ly + I2)/2 

Al = (Ul, — Iz)/2 

Tp = moment of inertia of blade about flapping 
hinge, Ibs. ft. sec.? 

Ir = product of inertia of blade, Ibs. ft. sec.? 

K, = spring restraint about flapping hinge, ft. lbs. 

Ky = spring restraint about feathering hinge, ft. Ibs. 

M, = aerodynamic moments about feathering axis, 


ft. Ibs. 


Ma 


Xb, Yo, 2b 
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aerodynamic moments about flapping hinge, 
ft. Ibs. 

inertia moments about feathering axis, ft. lbs. 

inertia moments about flapping hinge, ft. lbs, 

moments about flapping hinge due to spring 
restraint, ft. lbs. 

moments about feathering axis due to spring 
restraint, ft. lbs. 

moments transferred to helicopter by rotor, 
ft. Ibs. 

components of 1, and VM, due to inertia forces 
= 

moments produced by Ho to Hy 

H,(mRh/I) = H,(unyn/r) 

component of M, and M, due to Ky, = nww,?/2 

rotor radius, ft. 

real parts of Ay to A33 

amplitude of helicopter response to harmonic 
excitation” 

amplitude of servo response to harmonic exci- 
tation 

total disc area, ft. 

total rotor thrust, Ibs. 

tangential velocity at blade element, ft. sec.~! 

helicopter gross weight, lbs. 

coefficients of applied control moments defined 
in Eq. (20) 

slope of the lift curve of blade element, assumed 
to be 6 

blade chord, ft. 

radial distance from blade center of gravity to 
flapping hinge, ft. : 

flapping hinge offset from center of rotation, ft. 

Laplace transform 

acceleration of gravity, ft. sec.~! 

distance from rotor hub to helicopter center 
of gravity, ft. 

servo gear ratio 

characteristic time lag of automatic control 
system, sec. 

radius of gyration of blade divided by (R — e) 

helicopter mass, Ibs. ft.~! sec.? 

blade mass outboard of flapping hinge, Ibs. 
ft.—! sec.? 

Integral defined as dm, 

total number of blades per helicopter 

elements of roots of characteristic equation 

distance from center of rotor to blade element, 
ft., or servo gear ratio 

time, sec. 

velocity normal to blade element, ft. sec.~! 

disc loading, Ibs. ft.~? 

displacements of helicopter in X and Y direc- 
tions referred to axes fixed in space, ft. 

x+iy,g =x — iy 

displacement of blade element in X, Y, Z 
directions referred to axes fixed in space, ft. 

absolute angle of attack of blade element at 
azimuth position y, referred to axes fixed in 
space 

angles of helicopter pitch and roll, referred to 
axes fixed in space 

blade flapping angle at azimuth position y re- 
ferred to axes fixed in space 

coning angle of blades 

component of By in XZ plane 

component of By in YZ plane 

Bi + B = B, — 


a — B 
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nondimensional parameter defining rotor blade 


= 
characteristics = pabR*/I, = poaSR?/nIy, 

r = 

A = characteristic equation 

é = delta three effect such that 6 = tan 4; 

€ = nondimensional parameter defining helicopter 
characteristics = eopasR?/I = méy 

rd = nondimensional parameter defining rotor blade 
characteristics = emd/Ip 

= £+ + 6B 

n = nondimensional parameter defining rotor height 
above helicopter c.g. = h/R 

Oy = blade pitch at any azimuth position y 

6 = collective pitch setting 

61, 92 = cyclical pitch settings from trim position 

‘ = nondimensional parameter defining blade mo- 
ment of inertia = J,/I 

p = frequency variable of characteristic equation 

A, Az, As =.roots of characteristic equation 

Ae = rotor axial inflow divided by QR 

fe tte = ratio of helicopter forward speed in X¥ and Y 
direction to tip speed 

= 4/2 

~ = nondimensional parameter defining flapping 
hinge offset = e/R 

p = density of air = 0.002378 Ib. sec.? ft.~4 

o = rotor solidity, ratio of blade area to disc area 

T = nondimensional parameter defining helicopter 
characteristics = cpaSR/m 

$1 = phase shift between helicopter response and 
forcing furiction 

o = phase shift between servo response and forcing 
function 

— = chordwise distance between feathering axis 
and c.g. of blade element, ft. 

XA = chordwise distance between feathering axis and 
a.c. of blade element, ft. 

y = blade azimuth position at any time, ¢ 

Q = rotor speed, sec.~! 

Wb = natural frequency of restrained blade about 
flapping hinge divided by 2 

w = frequency of forcing function, sec.~! 


(II) THEORETICAL DEVELOPMENT 


The equations of motion will be derived by consider- 
ing the displacement of the helicopter and its blades 
relative to a system of axes fixed in space with the 
origin at the undisturbed trim position of the rotor 
hub (6; and ( initially zero). The equations will be 
linearized in the usual manner by assuming small dis- 
turbances from the equilibrium position. Second and 
higher order quantities may then be neglected, includ- 
ing such effects as a rise from blade motion in the hori- 
zontal plane. The inflow velocity will be assumed 
constant over the rotor disc and given by the customary 
expressions for thrust and torque in hovering flight: 


— 
(oC,,/8) + orCr 


Cr 
Co 


The theory will be developed in this section for the 
general case of an offset flapping hinge, but the results 
are directly applicable to a teetering rotor, when £ is set 
equal to zero, and to a rotor with blades rigidly at- 


Fic. 11. 


Coordinate system for rotor. 


tached to the hub, when the equations are modified as 
indicated in section ITI. 

Referring to Fig. 11, x is the displacement at any 
time ¢ of the rotor hub in the X direction and a, 6, the 
corresponding angular displacements of the helicopter 
and tip path plane, both measured from the trim posi- 
tions. Similarly, y, Be, and a» define the displacements 
in the Y direction. Angular displacements are positive 
in the direction indicated—that is, the thrust vector is 
tilted to the left (in the X direction) and out of the 
plane of the paper (in the Y direction). 

The equations of motion for the helicopter then be. 
summation of horizontal forces, 


come: 
m(% — + H, = 0 (la) 
— + H, = 0 (1b) 

summation of moments about the helicopter c.g., 
Ijin + Thea, + M, + hH, = 0 (1c) 
Thag hH, == 0 (1d) 


where m is the helicopter mass, including the blades. 
I, and J, are the moments of inertia of the helicopter 
(without blades, since moments due to the blade inertia 
effects are included in M, and M,) about its center of 
gravity. h is the vertical distance from the rotor hub 
to the center of gravity. ‘H, is the horizontal force in 
the X direction primarily due to rotor blade drag and 
the inclination of the thrust vector, including the in- 
clination resulting from a combination of coning and 
the variations in lift of individual blades with azimuth 
position. T is the total thrust, assumed constant and 
equal to the gross weight, since variations in thrust due 
to small order disturbances from the hovering condition 
appear only as functions of the products of the small 
order quantities 6, and yw, or and M, and M, are 
the moments about the X and Y axes exerted on the 
helicopter by the reactions at the flapping hinge when 
offset from the center of rotation and include the effects 
of both inertia and aerodynamic forces. Rotor head and 


| 
Yb 
- 
ec.—! 
lefined 
sumed 
diree- = 
Y, Z 
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fuselage parasite drag and fuselage damping moments 
are functions of second order quantities when small dis- 
turbances from the hovering condition are being con- 
sidered and therefore do not appear in these equations. 

The equations of motion for the blade may be de- 
veloped as follows: The position at any time ¢ relative 
to the fixed axes of an element of the blade of mass dm, 
and chord } at a radial distance r from the center of ro- 
tation will be (Fig. 11) to the first order of small quan- 
tities, and since y = M1, 


2 = (r — e)By + ea, (2) 
= x — rcos (3) 
=y—rsin Ot (4) 


JOURNAL OF THE AERONAUTICAL SCIENCES—AUGUST, 1948 


The moment about the flapping hinge due to the 
inertia forces, positive when it tends to decrease By, is 


M, = + (% cos Ot + Jp sin Mt)B,](r — e)dm, (5) 
The blade flapping may be expressed as 
By = Bo + Bi cos Qt + Be sin Mt 


higher harmonics of flapping having no effect on the 
stability of the helicopter as a whole. 8, and 2 are 
functions of time. {» is constant since the thrust is 
constant. 

Similarly, 


ay, = a, cos Qt + a sin MN 


where a; and a» are functions of time. 


Differentiating Eqs. (2), (3), and (4) twice with respect to time to obtain accelerations, substituting in Eq. (5), 
and performing the integrations, the following expression is obtained for the moment about the flapping hinge due 


to inertia forces, positive when it decreases B,, 


where J, = 


+ — 208; + Bs + — — + of a | sin + 


moment of inertia of the blade about the flapping hinge— 


I= L*(r — e)*dm 


and ¢ = distance of blade centroid from the flapping hinge divided by the square of the radius of gyration about 
the flapping hinge and multiplied by the hinge eccentricity e— 


= (e/h) — e)dm, 
The moment about the flapping hinge due to an elastic restraint of spring constant K, ft.Ibs. per rad. is giver 


by 


(M,/TIp) = 2? Bo + wp? (Bi a) cos Qt w»?( Bo a2) sin 


(7) 


where w, = VK »>/1p/Q, the natural frequency of the restrained blade about the flapping hinge, when 2 = 0, divided 


by 2. 


The moment about the flapping hinge due to aero- 
dynamic forces acting on the blade, positive when it 
opposes the motion, is 


M, = + uU)(r — e)dr (8) 


where U and u are velocities relative to a blade element 

dr, at radius r, and normal to the blade. Uisina plane 

parallel to the XY plane; w is in a plane perpendicular 

to the X Y plane and passing through the blade radius. 

Referring again to Fig. 11, 

—%, cos By — (%, cos Qt + % sin Qt) sin By, — 
ARQR cos By 


Substituting from Eqs. (2), (3), and (4), and since py 
is a small angle, > 


4#= (r — e)B, e&y — cos Ot — 
JB, sin QE — AOR (9) 


ARQR is the inflow into the rotor, assumed constant and, 
to the first order of small quantities, is the same whether 
considered normal to the rotor disc or the reference 
plane. 


U = sin Qt — Qt = X sin Qt — pcos Qt + Or (10) 


The blade pitch 6, is given, relative to the reference 
plane, as 


by = (A 5 Bo) + (0; sin Qt (05 + Q2) > 
cos — — a) cos Qt — — az) sin Qt (11) 


where 6) is the collective pitch and 6,, 6 the applied 
cyclic pitch variation. £8; and f are of the order of 


6, and 6, which establishes the order of 6, and 6, for 
which the equations are valid. 
hinge such that 6 = tan 43. 
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Substituting Eqs. (9), (10), and (11) im Eq. (8) and performing the integration, neglecting the effects of higher 
harmonics and higher order quantities as before and dividing by J,, gives 


M Dt C 


The coefficients are represented by the following integrals, which have been evaluated for the case of a constant- 
chord blade, and no tip loss factor: 


D= ab[2(% — 5Bo)(r — e)rR — dApR*(r — e)]dr = HEG — 65Bo) (5 


The equation of equilibrium for moments about the flapping hinge is 


This equation will be satisfied at any time ¢ if the sum of the coefficients of sin Qf, cos Qt and unity are sepa- 
rately equal to zero. 

Equating these coefficients to zero yields two equations describing the blade motion and an independent equa- 
tion from which §) may be determined: 


“Dap t + Ch Aa + + CES — (wy + — 
Q 


+ (wy? + + + + Ab, = + Bb (13b) 
Boll + £] + = — 580) — Crs (13¢) 


The forces and moments appearing in Eqs. (la) to (1d) may be developed as follows: 
The average moment M, about the Y axis due to restraint and offset of the flapping hinge, positive when it 
opposes the motion, is 


2x R 
M, = {1 + (% cos + sin Qt)a,|dm, — cos Qtd(Qt) + 


n 2x 
— — B) cos d(Qt) 


from which 


+ 205 + (or — + [ + ~ + Pe= + Py 


I 
mat 
— 68; + 6a) + (a a» | + 7 la 7 Loa & (14) 
Similarly, 
= 7 Bi + (a2 + | Ato + a + + OR 
Bo + + (ae a») TlLa = + T — & (15) 


M,+ M,+M, =0 (13) 

= 
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In these expressions the following integrals appear and have been evaluated for’ the case of a constant chord 


blade: 
L = (ep/4I) fc*abr* dr = (€/12)(1 — &) 
N = (ep/4l) ab(r — e)r dr = 
P= (ep/41) abr dr = (e/8)(1 — &) 
Q ab[2(8 — 5&)rR — dR? 


where = epoaSR?/I and I = (I, + J,)/2. 


The average horizontal force in the reference plane in the X direction, positive when it opposes the motion, is 


— pb 
me -if P 


and in the Y direction 
2a 
=p 


from which the H forces may be derived in the form: 


m 
Hy, + Hy (14a) 
& 


No inertia terms appear in these equations since, in 
the case of the horizontal forces, the blade mass may be 
included directly in the helicopter mass in Eqs. (la) 
and (1b). 

The integrals represented by the coefficients Hy to Hy 
have been evaluated for the case of a constant chord 
blade and are given below: 


rf 2 W T 

= 4L3 — — — m22R mR 

Bo? 
Hs = (r/4)[(QA/2) — §)] = — M, when 
He = (7/4)[(80/3) — 
Hy = (1/4) [(Bo/2)(1 — &) 
He = (1/4) — &) — 3A(1 — 


Hy 


(r/4)[1/s — (&/2) + 
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— 88)(1 — — 8)] 


An examination of these expressions for the forces 
and moments acting on the blades and the helicopter 
indicates that, if the variables * and 7 are replaced by 
their more usual forms uy, = */QR and p, = ¥/OR and 
if the differential equations of motion are reduced to 
algebraic equations by means of the harmonic substi- 
tution uw, = uz,e"", a = ae etc., then the coef- 
ficients of the characteristic equation will contain only 
the nondimensional parameters 7, 7, €, and dg, 
%, ¥ 

The labor involved in expanding the six-row deter- 
minant obtained from Eqs. (la) to (1d) and the two 
equations for the blade motion Eqs. (13a) and (13b) 
may be considerably simplified if the angular and linear 


displacements are replaced by their gay 8 forms: 
z=x+ 
a=at 1a, a= — 
B=Bitif, B= — 


and the moment of inertia in pitch, J,, and in roll, J, 

by J = (J, + J,)/2 and AI = (J, — I,)/2. This 

method of approach is suggested by Coleman's treat- 

ment of helicopter ground resonance in reference 1. 
The characteristic equation is then given by 


AA + (+, 4 — — = 
(16) 
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where 
Ay Ay» Ais 
A= |An Ax Ags 
Az Ax. Ags 


A is the conjugate of A and may be obtained directly 
from the expansion of A by changing the signs of the 
imaginary parts of the coefficients of the different 
powers of y. The coefficients Ay to As; are, neglecting 
the last two terms of Eqs. (14) and (15), which are 
negligible when « is small, as is always the case, 


Ay = (v+ — 

Ay = + Hw — (Ao — — i[Hw — Ap) 

An = [PB + M4] — z[(Q + Ms] 

Ay = [v? + o(PE + M7) + (MM + + My — 
— i[vms — (N + M)] 

Ax = + + Mo) — (M + + — 
Lé)] — i[v(2M + Me) + (N + 

An = {Bol(s/é)» + i 

An = [fot + Civ — + + 8B)] — i(2tv — A) 

Azz = [v? + Av + (€ + w? + 5B)] — i(2v + A) 


and A, to Ag; are again obtained by changing the signs 
of the imaginary parts of Ay, to Ags. SM; to My are ob- 
tained from H, to Hy from the relationship %, = 
H,(nym/r). Also Wg = and Mis mf/2. The 
nondimensional coefficient ¢ is based on I and not J,. 

The expansion of the six-row stability determinant 
then consists essentially of the expansion of a third 
order determinant and multiplication of the resultant 
fifth order equation by its conjugate form, giving a 
tenth order equation in v with real coefficients. It will 
be found that this equation is not readily solved di- 
rectly by ordinary methods, since the roots are all com- 
plex and appear in pairs of the same order of magnitude. 
An indication of the degree of stability, or instability, 
and the frequency of these oscillations may be obtained, 
without solving for the roots, from the frequency re- 
sponse and phase shift as discussed in section III; how- 
ever, the manipulation of the tenth order polynomial 
iseven then a tedious process. It will be the purpose of 
the next section to simplify the theory as much as ap- 
pears consistent with engineering accuracy. 


(III) APPLICATIONS 


Stability Characteristics 


Dual Rotor Helicopter —The first case to be con- 
sidered will be that of a twin rotor helicopter with 
either side by side, or coaxial, opposite rotating rotors. 
For this case the ‘equations of motion are reduced to 
four, since the lateral inclinations of the two tip path 
planes are equal and opposite and the rolling and pitch- 
ing oscillations of the helicopter are consequently de- 
coupled. 

Referring to section II, the equations of motion for 
this case, after converting to algebraic form by the 
harmonic substitution, uz = ure", a = ae’, etc., 
become: 


Ruz, + Ryan, + Rishi, + T1382, Ash, + 


Roz, + Raa, + + = + 

(9% + L)62) (17) 
+ Reo, + + = +Bé2 
+ T3201, + T3381, R3382, = —Bé, 


where the R’s and /’s are the real and imaginary parts 
of the coefficients Ay, to A33 listed in section II. Thus, 


Aun = Ru ily 


etc. The nondimensional coefficients « and ¢ for this 
case are based on J, and not J. 
For a helicopter with the following parameters, 


= 0.184 % = 0.159 A=0.050 =0 6=0 
y= 12 « = 0.029 «= 0.069 =0.5 n=6 
0.033 ¢ = 0.066 Bo = 0.14 = 0.012 


and 2 = 25, the roots of the characteristic equation 
obtained by setting the determinant A, represented by 
the left-hand side of the above 5 sem of equations, 
to zero are 


—18.35 + 8.311, 
—18.8 + 42.54 


A = —1.141, +0.2325 = 0.6771, 


It will generally be found convenient to convert the 
characteristic equation to time uriits expressed in sec- 
onds by the relationship \ = v@ before extracting the 
roots. 

The last two pairs of conjugate complex roots are pre- 
dominantly blade motions and represent heavily 
damped high-frequency oscillations, with a reduction to 
one-half amplitude in 0.04 sec. Consequently, it is 
evident that blade response following a disturbance 
will be rapid and will not be the cause of any delay in 
the response of the helicopter. The first pair of com- 
plex conjugate roots represent the unstable oscillation 
of the helicopter with a period of 9.3 sec. doubling its 
amplitude in 3 sec. and, together with the real root, 
establishes the response of the helicopter. 

The first step in simplifying the equations of motion 
will be an attempt to approximate the blade motions 
by the simple convergence represented by the real part 
of the last two pairs of roots, since their order of mag- 
nitude is such as to suggest that their influence on the 
helicopter motions will be small. This may be done 
by neglecting all direct acceleration terms proportional 
to ,, and also & in the equations of motion for the 
blade. Also, since £2 is small compared to 8; and the 
unstable root is small, 8; may also be neglected. With 
these assumptions the following coefficients will change 
to 


Is = Hy 

Roy = + Me) — (M + + My — Ld) 
= N+ M 

Rw = Civ 

R33 = Av + (ow? + £ + 6B) in the third equation 
Is, = A in third equation only 

Rss = (w? + ¢ + 6B) in fourth equation only 


The rest of the coefficients remain unchanged, 
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The expansion of the corresponding determinant gives a fourth order characteristic equation of the form 
Cw + Cy? + Cw? + Cv + Cs = 0 (18) 
whose coefficients may be approximated to a high degree of accuracy by 


CG = All + 
0 
= li + ({) Ke + P + a) | + 
[on +o + + te + (1+ Se) 
= 0 
= + + + 9% + am) + ay | + + (4) - + 


These coefficients may be further simplified when ¢ is small and w? = 6B = 0 to 
2, CG; =A, G=0 


where § = ¢ + 5B + a,?. 


D 
Cy = (ot + 914)(2 + wD) +4 4h] Pp, + — (N + | + 


Cs =, + + Ps + — (N + Me) 


For the case considered above this form of the sta- Bol 1 + . 2 _ Bzo ¢D eres 
bility equation gives the roots 20 i | ba 0 


=1.107, + 0.191 + 0.606, — 17.375 (1 
A 


which are in sufficiently close agreement with the com- 
plete solution for practical application. 

This is the simplest form of the stability equation 
which does not sacrifice generality. However, in the + (19) 
case of a conventional helicopter with small or no servo 
gear ratios such that the natural tendency of the tip 
path plane to follow pitching of the helicopter is not a 
restrained, further simplification to a third order equa- nities dhorcallewneianalial 
tion, similar to Hohenemser’s treatment in reference 2, + + a+eot+ (1 + + 
is possible by eliminating completely the rotor degree A[l + (¢/A)?] A 


The corresponding characteristic equation may be 
closely approximated by 


where 


of freedom and substituting for it a “quasi static’’ solu- (N + Me) [ 2- CE 
tion based on the assumption that, since differs (i +f) — (1 + 
from a by only a small amount Af; = a — fi, 6; may 
be assumed equal to The equations of motion + Mx 
then become, after reduction to algebraic form by the A[l + (¢/A)?] A 
substitution as before, (N + 91s) 
9 

= — — BC 
Ru uz, + (Ra + Ris)o1, + Hi + = 0 A[{l+ (¢/A)*]_ A 
X When is small and 6 = w,” = 0 these simplify to 


AB, + (N + = 0 
where K = (2/A)(M + 9%) 


1 + + J= (m+ om) + M+ a P gin + am) | 
Vai, cé\. For the case considered above this form of the sta- 
+0 + pitity equation gives the roots 


and A = —1.0, +0.185 + 0.67 
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which is evidently a close approximation to the de- 
sired roots. The stability determinant of reference 2 
expressed in the form used here becomes 


i 2 
A 
nyCp,m 
4 
whose roots are for the case being considered 


A = —1.125, +0.155 + 0.57% 


vt 
= 0 


2 


9% is derived from Hp, as were all other values of M,,. 


nm D D 


if the horizontal force along the tip path plane, 
(7/4) — (60/3)] 


due to “‘induced’’ drag, and > effects are neglected, as 
was done in reference 2. The main difference, apart 
from these small effects, between the two solutions is 
therefore the appearance of 91%; instead of S% in the ex- 


to the tip path plane, then the rearward inclination of 
the tip path plane Af; = (2/A)d, due to flapping 
arising from the Coriolis acceleration induted by the 
pitching velocity d:, would give rise to the damping 
moment 

Kéy = (2/A)(M + 


where IN is the inertia moment due to the eccentricity 
of the flapping hinge. However, the individual lift 
pf each blade element acts perpendicular to the wind 
direction relative to that blade element and, hence, is 
perpendicular neither to the tip path plane nor the 
reference plane. Consequently, feathering in the tip 
path plane or its equivalent flapping in the reference 
plane will always give rise to a horizontal force in 
either plane after integrating and averaging the com- 
ponents of the elemental blade lifts, in either plane, 
over the disc area. This has the effect of reducing the 
damping moment to 


= (2/A)(% + NM) 


Since this is a relatively small effect, it may be expected 
that, when £ is small and the rotor system otherwise 
conventional, the approximate solution given in refer- 
ence 2 will be in close agreement with the theory de- 
veloped here. The discrepancy between the results 
given in references 2 and 3 arises from the assumption 
in reference 3 that the blade pitch setting % and the in- 
flow ratio Ag are small order quantities, which has the 
effect of eliminating the rearward tilt, D/A, of the tip 
path plane due to forward motion. 

Single Rotor Helicopter.—Considering next the case 
of the single rotor helicopter, it would appear from the 
above analysis that a close approximation could again 


If the thrust were always perpendicular. 
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be obtained by making the same assumptions used in 
reducing the characteristic equation to the form given 
by Eq. (19). 

Using the complex form of the stability equation 
given in section II, the characteristic equation may be 
closely approximated by 


AA — (Al/I)*%* = 0 
where 
A= [(K + J) — (E + + — iF) 
and 
A= + [((K + oJ) + (E + + + iF) 
K and J have the same values as before and : 
E = M2 + Ms + 2M — 


[1 + L 
W+N 
All + 
m+ + De 
My + N [> 
Since rotors with large values of — would probably be 
used in pairs to minimize vibration, E and F may be 


simplified for the majority of present-day single rotor 
helicopters by assuming ~ to be small. Then, 


E = Mm — mM — (N + M,)(2/A) 


gt ato]- 


aa + + (1+ 


F=%+Q+ 


F=%+Q-—- (om + — (amy + 

In calculating K, J, E, and F for the single rotor 
helicopter the nondimensional coefficients « and ¢ must 
be determined using the average moment of inertia J 
and not J,, as in the case of the dual rotor helicopter. 

Table 1 lists the roots of the stability equations for 
a single rotor helicopter with various ratios of pitching 
to rolling moments of inertia and with the same non- 
dimensional parameters as were used in the case of the 
dual rotor helicopter. Such a helicopter would have a 
single rotor identical with one rotor of the dual rotor 
helicopter, half the gross weight, half the moment of 
inertia, and the same rotor height. If the single rotor 
helicopter were designed for the same gross weight as 
the dual rotor helicopter, the comparison would be 
slightly more favorable to the single rotor because of 
the relative increase in the nondimensional blade 
moment of inertia, 

It may be concluded from Table 1 that the single 
rotor helicopter will have a response in pitch somewhat 
similar to that of the dual rotor helicopter, since the 
unstable root is essentially unaffected by the rolling 
moment of inertia, as may be expected with such a 


arc] 
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loosely coupled system. However, the oscillations are 
no longer in one plane, and the appearance of a large 
roll component may be expected to further complicate 
the handling problem. If the helicopter were stabilized 
by means of a servo control to the extent indicated in 
Fig. 7, it should be possible to reduce these oscillations 
in roll following forward stick deflection to a small com- 
ponent by careful orientation of the control column 
and swash plate with respect to azimuth. However, 
in the case of an unstabilized single rotor helicopter, a 
large excitation in roll would appear to be inevitably 
associated with excitation in pitch, since the aerody- 
namic damping of the blade will always induce a tip 
path displacement at 90° to the initial displacement. 

Rigid Rotor—The equations of motion for a helicopter 
with infinitely rigid rotor blades fixed at the hub may 
be obtained directly from the equations of motion 
with an articulated blade by setting 6: = au, bi = du, 
and = 0 in the expressions for M/, and H,, substituting 
for M, the average aerodynamic moments on the 
blades, (1/21) f.°"M, cos y dy, and replacing J, by 
I,{1*+ (m/2)] since the rotors contribute half their 
polar moment of inertia to the moment of inertia of the 
helicopter. Since rigid rotors would only be used in 
pairs, Bp = a, = 0. If automatic feathering 6; propor- 
tional to w is incorporated in order to relieve blade 
bending moments at high forward speeds, Hi, will be 
decreased slightly because of the absence of any ap- 
preciable “induced drag’ in the tip path plane. It will 
be assumed here that no such feathering exists. 

With these modifications the complete stability de- 
terminant becomes 


My, 


—nv? + Hw — 


Using the same nondimensional parameters as be- 
fore, the corresponding characteristic equation has the 
roots 

A = —3.76, +0.032 + 0.567 


which shows only a slight improvement over the case of 
the articulated blade. 

Articulated Blade with Offset Flapping Hinges.— 
In order to establish the effects of a radical change in 
rotor design involving a large offset flapping hinge, a 
twin rotor helicopter similar to the one used in the 
previous examples but with £ increased ten times will be 
considered. It will be assumed that the blade moment 
of inertia about the flapping hinge is proportional to 
(1 — &)* and that the radial c.g. and radius of gyration 
are the same percentage of the flapping length of the 
blade. Since the fixed portion of the blade will con- 
tribute practically nothing to the rotor lift, it will be 
assumed that the aerodynamic characteristics of the 
rotor remain the same. The following nondimensional 


parameters will then change, the rest remaining the 
same: 
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= 0.33, y= 40, ¢ = 1.0, « = 0.00875, = 0.69 


With these values the roots of the stability equation 
obtained from the complete determinant for a twin 
rotor helicopter become 


A = —6.55, +0.0294 += 0.5761, —18.1 = 1.87%, 


—20.9 + 53.54 
Eq. (18), the simplified form of the stability equa- 
tion, gives the roots 
= —5.76, +0.04 +0.58i, — 14.72 
and Eq. (19), the most simplified version, gives the 
roots 
A = —4.287, +0.033 = 0.582 


which is adequate for most design purposes and may be 
uised for all cases except when the helicopter and rotor 
are decoupled in pitch. 


Control Characteristics 


The control response of the twin rotor helicopter 
only will be considered, since the stability equations 


indicate that the single rotor helicopter may be ex- 
* pected to have similar control characteristics in any one 


plane. 

The system of equations represented by Eq. (17) 
may be solved for any of the degrees of freedom, a, 
41, 61, B2 by the usual methods of operational calculus. 
If, instead of the harmonic substitution e’™, the opera- 
tion f,°e~’™ f(t)dt is used in solving the differential 
equations of motion, then the Laplace transform L£(a,) 
for the pitching degree of freedom following unit step 
function applied to 6; and 4 is given by ‘ 


_ + ( 1 ) 


(an) A vQ 

Since it will be found, in general, that the large con- 
jugate complex roots representing blade response do not 
affect the response in pitch to any noticeable extent, A 


may be used in the form given by Eq. (19). The cor- 
responding values of G(v) is G(v) = Xv + Z, where 
Afi + A*[1 + 
and 
BH, 
Z=- = — 
A[i + (/A)*] 
Also g(v) = Yv + V, where 
BIN + Ms) , FB(M + 
Y= -— = +L+ 
A[l + (§/A)*] A®[1 + (¢/A)?] 
and 
BH, (BH, 
V=- = + H, 
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When é is small, these simplify to 


X = —(M + M + Ms) = —(M + Md) 
Y=0 


In determining initial response, Z and V may always be 
neglected. 

The inverse transform of £(a,) may be obtained di- 
rectly, by expanding in partial fraction, in the form 


e™ 


— Ar)(As — Aa) 


= (0X + 


+ | 
(Ar — Ae)(As — Az) — As) (A2 — As) 
where Ai, Az, As are the roots of the characteristic equa- 
tion 
A = (A + + Az)(A + Az) 

If these roots are written A; = s, A» = p + ig, As = 
p — iq then the transient response of the helicopter be- 
comes 

a 


= — 
(Xi + YR) +q? 
cos gt — : sin a | (20) 
For the dual rotor helicopter with small flapping 


hinge offset — = 0.033, Eq. (19) gives, in second units, 
s= 1.0, p = —0.185, gq = 0.6. Also, 


OX = —0.0205 2? = —12.8 


This case has been plotted in Fig. 4 and discussed in 
section I. 

A more accurate representation of the transient re- 
sponse can be obtained by using the pertinent roots 
obtained from the complete solution, when available, or 


s=1.141, p = —0.2325, q = 0.677 


Fig. 9 shows the response when & is large (£ = 0.33). 
For this case, Y is not zero and contributes appreciably 
to the control moment. If the control column linkage 
were arranged in such a manner that forward stick de- 
flection corresponded to initial forward inclination only 
of the tip path plane with no lateral inclination, then 
= —(¢/A)@, or, more generally, 


— + w? + 


since, in the static condition, 


Ey) _B 
E + ( A B, = A 2 + A 1 
For the case plotted in Fig. 9 this condition is satis- 


fied when 6. = —0.66;. The values of X and Y used 
in Fig. 9 are, in second units, 


= 


—46.5 
29.4 


= 0.07440? 
2?Y = 0.0472? 
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SERVO GEAR RATIO, k 


—— Approximate Theory Based on Equation!8 
—-— Approximate Theory Based’on Equation 
—-—--Exact Theory Based on Equation !7 


Fic. 12. Effect of servo gear ratio, k, on oscillatory mode of 
helicopter. ‘ 


For the case of the rigid rotor plotted in Fig. 8, Eq. 
(20) for the motion following unit control deflection may 
be used if X is set equal to 0 and Y is replaced by 


= 


in second units, which, for the case being considered, 
gives 2?Y = 82.5 or over six times the control power of 
the flapping rotor. The comparison is not, however, a 
complete indication of the relative control effective- 
ness, since, as discussed in section I, a considerably 
higher cyclically applied blade deflection could be 
tolerated in the case of the flapping rotor. 

Control Response with Automatic Stabilization.—li a 
gyro stabilizing unit is incorporated in the helicopter 
control system in such a way that a cyclical pitch con- 
trol 4, is applied to the rotor so that 6; = k(a: + ra), 
where & and kr are the servo gear ratios, then Eq. (18) 
takes the form 


Ct + [OC, + + + + 
A(M + My) + [A (MM + IM) + = O 


and Eq. (19) takes the form 


+ [Q(K + nJ) — — [Q°RX]A + = 0 
(21) 


It should be emphasized that Eq. (21) is only good for 
small values of k. Fig. 12 shows the variation of the 
real and imaginary parts of the complex conjugate root 
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with k, when r = 0, as determined from the complete 
solution, and the approximate solutions corresponding 
to Eqs. (18) and (19). It is evident that the maximum 
degree of stabilization is obtained using a small value 
of k, and, consequently, it is believed that Eq. (21) is 
adequate for most applications. In establishing the 
transient response, Eq. (20) is again directly applicable 
when. the values of s, p, and g obtained from Eq. (21) 
are used. The transient response for various combina- 
tions of r, k, and & have been plotted in Figs. 6 to 10 
and discussed in section I. 

In practice, the displacement signal, ka, may be ob- 
tained by means of a displacement gyro or by electrical 
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integration of the rate signal from a restrained gyro. 
The latter arrangement has several advantages but 
would have to incorporate a ““memory”’ if its reference 
is to remain unchanged with either helicopter or control 
surface deflection, as has been assumed in this analysis, 
Similarly, the rate signal (rk)a@ could be obtained 
either directly from a rate gyro or by electrical differen- 
tiation of the displacement signal, ka. In either case 
time lags are involved in both the servo system and, to 
a lesser extent, in the restrained gyro. In conventional 
aircraft autopilots, these time lags have been of the 
order of 0.1 sec., although by careful design this figure 
may be reduced appreciably. 


If a characteristic time lag / exists in the servo system such that 


16; + A; = 


then Eq. (18) becomes, for small values of §, 


k(a + ra) 


+ [Ci + + + 10°C; + 2MskrO?]A¥ + + + A(M + + + 


+ + = 0 


and Eq. (19), applicable only when & and / are small, becomes 
+ [1 + Q(K + + [Q(K + nJ) — — — + = 0 (22) 


For any but the simplest servo system it will usually be found preferable to treat the characteristics of each ele- 
ment separately and combine them by the frequency response method discussed later. 
Flight experience has shown that a considerable change in helicopter control characteristics may be realized 


by shifting the blade chordwise center of gravity location relative to the feathering hinge. 
change may be established by an examination of Eqs. (6) and (12). 


The effects of sucha 
If the integrations (5) and (8) are performed 


with (r — e) replaced by the distance x, between the feathering hinge axis and the chordwise c.g. and the distance 
x, between the aerodynamic center and the feathering hinge axis, both assumed positive for displacements ahead 
of the feathering axis, then, neglecting the effects of Coriolis accelerations in the plane of rotation, the inertia and 
aerodynamic moments M, and M, about the feathering hinge would be, positive when tending to decrease @, 


= Si" + (a, cos Ot + jp sin 24) By]x,dm, 
M, = —(*/2)p ab(U%, + uU)x,dr 


where f is the distance between the center of rotation and feathering hinge. 
the feathering axis is assumed to be negligible. 


The blade moment of inertia about 


The equations of motion about the feathering hinge C’ = (pab/2Ip) J? Rrx,dr 
may thus be derived, in the same manner as Eq. (13), as B’ = (pab/2Tp) Six dr 
M,+ M,+ M, =0 = pab/2Ip — 68) — 


where M, represents the control restraint about the If 
feathering hinge. M,and M/, are similar in form to M, 
and M, derived in section II except that in place of J, 
the product of inertia, 


M, = Kp + 


then, if a2, u, are zero, 


Tr 2 Tp 
would appear and, in place of /,¢, 
I, = eS;*x,dm A (a1 — Br) + D'u, + BY — be (23) 
and. 2. would, be where 6; is the sine component of the resultant cyclic 
modified to 


pitch variation and the damping force cé is applied to 


the nonrotating element of the swash plate. Positive 
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values of x, thus correspond to an increase, and nega- 
tive values to a decrease, in rate damping. Conse- 
quently placing the blade c.g. ahead of the feathering 
axis may be expected to improve the helicopter sta- 
bility characteristics, at the expense of control effec- 
tiveness, while placing the c.g. aft of the feathering 
axis will result in a quicker response at the expense of 
both blade and helicopter stability. Either case may 
be approximately simulated in Eq. (21) by an equiva- 
lent rk term, 2P'/Q, where = 

When TP is large, such an approximation is no longer 
valid, and Eqs. (17) or (18) should be used. Allowing 
for the effects of x,, Eq. (18) takes the form, for small 
values of f, and x, = J, = 0: 


CUS + M[C, + MC, + AT] + + + 
+ + Mo)] + mc, 


D 
AT(MN + am) ( A 1.) + C;24 = 0 (24) 
where the time lag / = c/K, arises from the swash 
plate dampers with damping constant c. When / = 
0 and L is 10, this equation, for the twin rotor helicopter 
with small hinge offset, has the roots 


A = 0.585 + 3.537, —0.009 + 0.1977 


The complete solution, as given by Eq. (17) and in- 
cluding both 4, and @, effects has, for the same condi- 
tions, the roots 
—1.13, —0.0073 += 0.19711, —0.635 + 3.571, 

—1.165 + 507 


1 = 


The effect of x, is thus to increase the period and damp- 
ing of the previously unstable root to a point where its 
effects on the initial control response would be negligi- 
ble. The blade response is no longer so rapid as before 
and becomes one of the dominant factors establishing 
the initial control response of the helicopter. This 
response will be poor, a/@, being of the order of 0.75 
when ¢ = 1 sec., and it would appear desirable to pro- 
vide some lag, /, by means of swash-plate dampers in 
order to avoid any appreciable reduction in blade re- 
sponse to cyclical pitch control. This will have the 
additional effect of providing some displacement re- 
sponse, by virtue of the integrating characteristics of 
the dampers, and also will damp out any oscillation due 
to blade inertia about the feathering hinge. If / = 1, 
and lr = 10 as before, Eq. (24) has the roots 


A = —26.96, —0.52 = 1.941, —0.0386 + 0.2451 


The vatue of a;/6,, neglecting the effect of the large 
root representing blade response and assuming the 
small roots to be zero, would then be a/@; = 1.5 at 
t = 1 sec., which may be further improved by increas- 
ing & providing J, is kept zero. The response curve 
would be similar to Fig. 2, and,.as.in the case of the 
phugoid oscillation of the airplane,,the helicopter would 


eventually return to the initial position by virtue of the 
long-period oscillation. 

The provision of a sufficiently large amount of T to 
realize such control characteristics involves several 
problems, probably the most serious being the twist 
of the blade in the steady-state condition due to a com- 
bination of coning, centrifugal force, and the offset 
gravity axis. This problem may be overcome with 
torsionally stiff blades by interconnecting the blades 
through a conventional swash plate with an irreversible 
collective pitch control. However, it may usually be 
advisable to balance the moments out directly at the 
blade by means of a combination of positive x4 and 
some camber in the blade. The equations of motion 
about the feathering hinge, Eq. (23), combined with 
Eq. (13a), then yield the following expression for the 
rearward inclination of the tip path plane relative to the 
helicopter shaft axis, with the simplifying assumptions 
that 6, = & = 8, = 0 and £ is small, 


D 


ab*R’ 1 A’ 

from which it may be seen that an increase in damping 
may be obtained by means of negative x, as well as 
positive x,. 

The summation of steady-state moments about the 
feathering hinge gives the equation 


AA 
Bo o — + 61, 


r 


where A is the change in collective pitch setting of 
the blade due to x,, x4, and C,,, the pitching moment 
coefficient of the cambered blade section. 

Evidently, if A®@ is to be zero, 


Bo = — + 


One satisfactory combination of x, x4, and C,, 
would be that which satisfies this equation and also 
reduces the rearward inclination of the tip path plane 
to 


which may be further modified to eliminate completely 
or even change the sign of the rearward inclination of 
the tip path plane due to forward. velocity u,, thereby 
eliminating the long-period oscillation. However, it 
is believed that this step is not advisable because of the 
possibility of divergence associated with such a modifi- 
cation and the desirability of eventually returning the 
helicopter to its initial position. 

With these values of x4, x, and C,,, Eq. (24) is 
again applicable when [ is multiplied by the factor 
[1 — (A’/A)]. 
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FREQUENCY- SECONDS 


Fic. 13. Helicopter frequency response and phase shift — = 
0.033. 


It is evident that the elastic characteristics of the 
blade in the chordwise direction may be expected to 
have an appreciable effect on the control response, par- 
ticularly with tubular spar blades, due to the shift in 
blade c.g. and aerodynamic axes, relative to the 
feathering hinge, associated with such deflections. 
Further second order effects may be expected from the 
Coriolis acceleration in the plane of rotation associated 
with flapping. A discussion of the elastic deflection 
in the plane of rotation and the forces producing these 
deflections is contained in reference 4. 

Frequency Response.—The frequency response 
method provides a useful means of establishing the ef- 
fects of various automatic control systems on the sta- 
bility of the complete loop made up of the helicopter, 
the various gyro elements, and the servo motors and 
has been described and developed in several texts and 
reports on servomechanisms. In particular, reference 
5 gives a simple introduction to this method with specific 
application to the aircraft control problem. 

The frequency response of the helicopter in pitch, 
Ry, taking for illustration the simplified form of the 
stability equation resulting in Eq. (19), is obtained 
directly by assuming a solution of the form e+” 
to the differential equations and solving for a/4. 
This gives 


iwXO? + (25) 


Z is, in general, negligible. 

For any specific forcing frequency, w, Eq. (23) will 
give a complex quantity, the absolute value of which is 
the desired amplitude ratio, or response. The phase 
angle ¢;, between the applied control deflection 6, and 


the helicopter response in pitch, a, is given by the ratio 
of the imaginary to the real part of this complex quan- 
tity. Fig. 13 shows the frequency response and phase 
shift for a helicopter with & = 0.033 and Fig. 14 that 
for a helicopter with = 0.33. 

The response obtained’ when Eqs. (18) and (17) 
are used to establish Ry are also shown in Figs. 13 and 
14 for purposes of comparison. The values of tan ¢, 
are plotted for simplicity in establishing the critical 
value of /, and only in the region where ¢ passes 
through 180°. The basic helicopter characteristics 
are the same as those used in the previous examples. 

Any desired servo response, obtained either experi- 
mentally or from theory, may be superimposed on Figs. 
13 and 14 in order to obtain a quantitative indication 
of the stability boundaries for the closed loop and a 
qualitative indication of the degree of stability and 
transient response. 

In the frequency range that is of primary interest for 
the helicopter, a conventional servo system may be ap- 
proximately synthesized by neglecting inertia effects 
and assuming the equations of motion to be given by 


+6 = + (26) 


when & and kr are the control gear ratios and / is the 
characteristic time lag corresponding to the servo 
damping factor, assumed viscous. Then the servo 
frequency response, kR,, is given by 
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Fic. 14. Helicopter frequency response and phase shift § = 0.33. 
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Combining Eqs. (25) and (27) gives 
kR,Rye — 1 = 0 


which is the same as A given by Eq. (22), previously 
developed by the method of determinants, when a solu- 
tion of the form e“” is assumed. This solution will be 
satisfied when A = 0 or 


wd 
which in turn will be satisfied if 
RR.Ry = 1 


when ¢2: + ¢: = 0. 

But a solution of the form e represents a boundary 
condition between stability and instability and, conse- 
quently, if, at any frequency w, ¢: + ¢: = Oand kR,Ry 
= 1, a condition of neutral stability will exist. 

If the servo system has no time lag (J = 0) and if 
r = 0, then ¢@ = 0 and R, = 1. Evidently, for the 
case of Fig. 13, a servo gear ratio k = 0.05, where Ry is 
taken at the value of w for which tan ¢; = 0, is sufficient 
to ensure a condition of neutral stability. This is in 
agreement with the conclusions of Fig. 12. 

If r = 0 but / is finite, then 


tan = —al 


and the value of / at which any amount of & will not 
stabilize the system occurs when the straight line repre- 
senting — tan ¢: becomes tangent to the curve corre- 
sponding to tan ¢:. When & = 0.033 (Fig. 13), this 
occurs at a value of / of 0.28 sec. from the approximate 
Solution or 0.18 sec. from the exact solution and when 

= 0.33 (Fig. 14) at a value of / = 2.7 sec. 

If r is finite, 


tan ¢d. = w(r — 1)/(1 + wl) 


and evidently, for ‘the order of frequencies encoun- 
tered with helicopters, small lags in the displacement 
control may be exactly compensated by rate control. ° 


APPENDIX 


Control Response of Conventional Airplane 


Pitch.—It will be assumed that the initial response 
in pitch is not appreciably affected by the long-period 
oscillation. With this approximation the equations of 
motion in pitch in nondimensional form® become 

7°76 + + aum, = bum; 
1A 
7(6 — — az, = 0 (1A) 


if the change in lift due to elevator deflection is ne- 
glected. The coefficients are 


(Ge) Ac) 
2kg\da/z\c/ S 


m, = (1/2k,)(dC,,/da) 


(=)! 


= 1/2(dC;/da) 


S = wing area, sq.ft. 

S, = tail area, sq.ft. 

¢ = mean aerodynamic chord of wing, ft. 

Ll = tail length, ft. 

m = airplane mass, lbs. sec.? ft.—! 

k, = B/mc? where B is the moment of inertia about 

the c.g. in pitch 

U = forward velocity, assumed constant, ft.—! sec. 

r = m/pSU, w= m/pSc 


6 is the airplane altitude relative to earth axes, and a is 
the angle of attack relative to the flight path. 


The transient response may be readily obtained by operational methods, as was done in the case of the heli- 


copter in section III, as 


= [um;/7?(p? + gt + sin (2A) 


where p and g are the real and imaginary parts to the roots of the characteristic equation obtained from the right- 


hand side of Eqs. (1A). 
Similarly, 


a —2p t 


(p* + ¢*)* 


( 2 cos gt + (sin | (3A) 


Figs. 1 and 2 have been plotted for a value of U of 146 ft. per sec. (100 m.p.h.), at sea level, and 


m, = 2.34, 


2 m, = 0.1, 


m; = 0.585, 


= 2.25 = 40, = 1.35 


Roll.—It will be assumed that the initial response in roll is not appreciably affected by the long-period 
convergence, or divergence, and the oscillatory root representing “Dutch roll.’ 


in- ; 
ise 
lat 
7) i 
nd 
cal 
ics 
gs. 
a 
‘or 
ts 
6) 
- wh 
he where | 
vo 
| 


472 JOURNAL OF THE AERONAUTICAL SCIENCES—AUGUST, 1948 


as Then : Fig. 3 has been plotted assuming values of 
= dL, (4A) L,=35, L,=120, U=100mph. 
where ¢ = angle of bank and ‘ Sectiiaiies 
Lp = C,,(gSb/I,)(b/2U) 1 Coleman, R. P., Theory of Self Excited Mechanical Oscilla- 
? tions of Hinged Rotor Blades, N.A.C.A. A.R.R. No. 3G29, July, 
Ly = C,,(9Sb/Iz) 


2? Hohenemser, K., Dynamic Stability of a Helicopter with 
Hinged Rotor Blades, Ingenieur-Archiv, Vol. 9, pp. 419-428, De- 
and J, is the moment of inertia in roll. f cember, 1938 (for English translation see N.A.C.A. T.M. No. 
If C,, = 0.5, C,, = 0.20, then the helical angle for 997, 1938). 

6, = +15° aileron deflection is given by 3 Donovan, A. F., and Goland, M., The Response of Helicopters 
with Articulated Rotors to Cyclic Blade Pitch Control, Journal of 

pb/2u = 0.10° the Aeronautical Sciences, Vol. 11, No. 4, October, 1944. 
* Yuan, Shao Wen, Bending of Rotor Blade in the Plane of Rota- 


a representative value; p is the rate of roll ¢, and bis tion, Journal of the Aeronautical Sciences, Vol. 14, No. 5, May, 
1947. 


the wing span. ibs Th 
From Eq. (4A) 5 Greenberg, H., Frequency-Response Method for Determination : 

ai of Dynamic Stability Characteristics of Airplanes with Automatic turbo 

Controls, N.A.C.A. T.N. No. 1229, March, 1947. blade 
eg — =x — | 6 Jones, B. M., Dynamics of the Airplane, in Durand’s Aero- of fini 
&, dynamic Theory, Vol. V, Division N. rotati 
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The Flow of a Perfect Fluid Through an 
Axial Turbomachine with Prescribed 


Blade Loading 


FRANK E. MARBLE* 
California Institute of Technology 


SUMMARY 


_The theory of the three-dimensional flow through an axial 
turbomachine, associated with variation of circulation along the 
blade length, is described as an extension of the classical theory 
of finite wings and is simplified to a problem in axially symmetric 
rotational fluid motion by considering an infinite number of 
blades in each row. The problem is linearized by considering 
the vorticity generated by the blades to be transported by the 
mean velocity. 

The linearized problem leads to well-known partial differential 
equations and is solved for the radial, tangential, and axial veloc- 
ity components induced by a single row of stationary or rotating 
blades with finite chord and prescribed loading. The particular 
case for which the blade chord approaches zero and the tangential 
velocity changes discontinuously is associated with the theory of 
the Prandtl lifting line for finite wings. Because the problem has 
been linearized, the solution for a multistage turbomachine 
follows by superposition of appropriate single rows. 

Analytical expressions and graphical values for the velocity 
components are given for a single stationary or rotating blade 
row of given loading with a hub/tip ratio of 0.6 and a blade aspect 
ratio of 2. The corresponding discontinuous approximation 
is compared with the more nearly exact solution and is shown 
to constitute a useful approximation to the solution for a 
finite blade chord when the discontinuity is located appro- 
priately. 

An exponential approximation for the velocity components, 
deduced from the analysis, allows rapid estimation of the rate at 
which the equilibrium velocity profiles develop ahead of, and be- 
hind, a blade row and, using the superposition principle, provides 
a simple means of approximating the velocity distribution in a 
multistage turbomachine and of discussing mutual interference 
of blade rows. 


INTRODUCTION 


a AERODYNAMIC PROBLEMS Of axial-flow turbo- 
machine theory which are subject to treatment by 
the theory of perfect fluids may be classified under the 
following two general problems: 

(1) Given the blade loading, blade speed, and the 


fluid state far ahead of all blades, determine the three- 


dimensional velocity field, blade shape, and distribu- 
tion of energy in the fluid. 

(2) Given the blade shape, blade speed, and the 
fluid state far ahead of all blades, determine the three- 
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dimensional velocity field, blade loading, and the dis- 
tribution of energy in the fluid. 

These are designated, respectively, the inverse and 
direct problems of turbomachine theory in analogy 
with the corresponding classical problems in the theory 
of finite wings. 

The inverse problem, which is definitely the less dif- 
ficult of the two, derives this advantage from the fact 
that it may be separated into two independent prob- 
lems. The first of these consists in determining the 
velocity field corresponding to the prescribed blade 
loading by replacing this blade loading with an equiva- 
lent force system distributed over some definite surface 
or region. In the second step, the two-dimensional 
theory of airfoils or airfoil lattices is applied to find 
the radial distribution of blade shape and orientation 
providing the prescribed load distribution when placed 
in this particular velocity field. Recent investigations 
into the theory of airfoil lattices (e.g., Garrick,’ Light- 
hill,? and A. Goldstein and Jerison*) have reduced 
most cases of the second step to one of numerical cal- 
culation. The present paper is concerned with a quan- 
titative description of the three-dimensional velocity 
field that prescribes the mean flow in which each blade 
element is situated and indicates the accuracy with 
which the two-dimensional airfoil theory may be applied 
to each element of the blade. 

For a perfect fluid the solution to the three-dimen- 
sional problem is simple only when the distribution of 
tangential velocity in any plane normal to the turbo- 
machine axis is that of a vortex situated on the axis.‘ 
Under this condition no radial or axial disturbances are 
induced, the circulation about each blade element is 
constant along its length, and the blade behaves sim- 
ilar to an infinite wing. In spite of the increased 
complexity of the flow under more general conditions 
of varying circulation along the blade, the approximate 
difference between the axial -velocity profiles far up- 
stream and far downstream of the blade row is easily 
calculated (e.g., Traupel,® Sinnette,* and Eckert and 
Korbacher’) by neglecting the radial transport of vor- 
ticity shed from the blades. From considerations of 
the simplified vortex system (e.g., Ruden®*), it follows 
that one-half of this axial velocity change has taken 
place by the time the fluid reaches the blade row. The 
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Fic. 1. Coordinate system and designation of velocity and 
vorticity components. 


Fic. 2. Deformation of a vortex ring. 


velocity distribution in the rest of the flow field, which 
is of importance in answering such practical questions 
as the interference of adjacent blade rows in a multi- 
stage turbomachine, cannot be found so simply. 

For the discussion of the three-dimensional flow in the 
present paper, the physical problem is simplified by 
considering a nonviscous and incompressible fluid, inner 
and outer boundaries consisting of concentric cylinders, 
and by assuming an infinite number of blades in each 
blade row so that the flow possesses axial symmetry. 
Therefore, the vorticity shed from each blade row is no 
longer concentrated in sheets but is continuously dis- 
tributed over the region downstream of the blade row. 
The difficulty of this problem lies in the nonlinear par- 
tial differential equations that describe the rotational 
fluid motion. By the approximation that the vorticity 
is transported by the mean velocity and not by its own 
induced velocity, the problem reduces to the solution 
of a well-known nonhomogeneous linear partial dif- 
ferential equation. The resulting solutions provide 
the linear approximation to the radial, tangential, and 
axial velocities associated with any loading of a blade 
of finite chord or of an infinitely thin blade row corre- 
sponding to a discontinuous change of tangential veloc- 


ity. 
FORMULATION OF THE MATHEMATICAL PROBLEM 


The flow is described (Fig. 1) in a cylindrical coordi- 
nate system 7, 3, z, by the velocity components u, 1, 
w, respectively. The corresponding radial, tangential, 
and axial vorticity components are 


= —0v/dz (1) 


n = (0u/0z) — (Ow/Or) (2) 


= (2/rdr) (or) (3) 


Because of the axial symmetry, only the tangential 
vorticity is associated with the radial and axial veloci- 
ties, while the radial and axial vorticity components 
are associated with the tangential velocity. Conse- 
quently, information concerning the radial and axial 
velocities is obtained by considering (cf., Meyer®) an 
annular vortex ring of small cross section (Fig. 2) which 
consists always of the same fluid elements. In par- 
ticular, the circulation F about the ring is given by I 
= 7 X cross-section area and remains constant through- 
out all deformations in a conservative force field. The 
only deformation of this ring consistent with the as- 
sumption of axial §ymmetry is a radial stretching, and, 
inasmuch as the volume of the ring must remain con- 
stant, the cross-sectional area varies as the ratio r/ro, 
where 7) and r are the initial and final radii. Under 
such a deformation the constancy of circulation de- 
mands that the initial and final vorticities, m) and 7, 
satisfy the relation 


n/r = o/ro = constant (4) 


Hence, as the ring moves along the stream surface in a 
conservative force field, the rate of change of n/r 
vanishes. 

The circulation about the cross section of the annular 
vortex will vary, however, in the presence of a non- 
conservative force field. Since the pressure is a scalar 
quantity and, consequently, generates a conservative 
force field, the circulation may be changed only by the 
action of the centrifugal force and the forces applied 
by the blades in the radial and axial directions The 
rate of change of circulation is easily found from the 
equations of motion: 


P 

= + 2(F) (5) 
Or \p 

wi—ut = —F; (6) 

un — vi = —F, (*) (7) 
Oz \p 


where P is the local stagnation pressure of the fluid and 
F,, Fs, and F, are the applied forces per unit mass in 
the radial, tangential, and axial direction, respectively. 
By subtracting the partial derivative by r of Eq. (7) 
from the partial derivative by z of Eq. (5) and simpli- 
fying by means of the continuity equation 


(Ou/Or) + (u/r) + (Ow/dz) = 0 (8) 

it follows that 
On un On OF, oF, 


Finally, by dividing by 7, simplifying the term on the 
left, and applying the definition of the vorticity com- 
ponents to the term on the right, it follows that 


where . 

The 
left sid 
tdken 
so that 
cylinde 
radial | 


If furt] 
small v 


then th 
tion to 
vorticit 
[Eq. (8 


or 


For 
in Eq. 
gential 
equatio 


in whic 
a strea: 
the tur 
tangent 
above, 


where | 
vorticit 
blade re 
velocity 


09(712) 


where 7 
of the fi 

The | 
lated tl 
act no! 


os 
’ 
\ 
| 


FLOW 
ds _ (u2 + 2)(2) 
dt \ o dz/\r 
1f oF, OF, 
9 
(9) 
where s is the distance measured along a stream surface. 
The operator u(0/Or) + w(0/0z) occurring on the 
left side of Eq. (9) represents the time rate of variation 
tdken while moving with the stream. If |u| << |w| 
so that the stream surfaces are appyoximately coaxial 


cylinders and if the tangéntial vorticity has a smooth 
radial distribution, then 


|u(On/dr)| << |w(dn/dz)| 


If furthermore the radial variation of axial, velocity is 
small with respect to the mean value w—i.e., 


| (w wo) / 1 


then this operator may be taken with good approxima- 
tion to be w(0/O0z). From definition of the tangential 
vorticity and by applying the equation of continuity 
[Eq. (8)], 


r dz 


r Or\r 02? 

For calculation of the centrifugal forces that occur 
in Eq. (9), it will be sufficient to approximate the tan- 


gential velocity in the following manner. The second 
equation of motion, Eq. (6), may be rewritten 


[u(0/Or) + w(0/Odz) = rFy (11) 


in which form it states that the rate of change, along 
a stream surface, of the moment of momentum-about 
the turbomachine axis is equal to the moment of the 
tangential force. To the same approximation as 
above, 


Fy ~ w(0v/0z) ~ wo(dv/O0z) (12) 
= —WEe 


where £ corresponds to the density of “bound”’ radial 
vorticity associated with a particular point of the 
blade row. From Eq. (12) the approximate tangential 
velocity becomes 


— v_.(r) = — J &(r,8)dB = 


F5(r,8) 
dp (13) 


where v_..(7) is the initial tangential velocity far ahead 
of the first blade row. 

The forces applied to the fluid by the blades are re- 
lated through the condition that their resultant must 
act normal to the blade surface. Furthermore, the 


OF A FLUID THROUGH AN 


AXIAL TURBOMACHINE 475, 


blade surface is parallel to the relative velocity of the 
fluid, and, consequently, the forces exerted by the 
blades are normal to the relative velocity of the fluid, 
If the blades move about the axis with an angular 
velocity w, this condition is expressed through the rela- 
tion 

uF, + (v — wr) Fs + wF, = 0 (14) 


But since the radial force is usually of smaller order 
than the tangential and axial forces, this becomes very: 
nearly 


(v wr) Fy + WF, 0 (15) 


inasmuch as 


Wo | | Wo | 


Now Eggs. (12), and (15) may be used to express the 
right-hand member of Eq. (9) in terms of known func- 
tions 


r Oz \r r \ Oz or) 
10 re) 
— of) (vor )+ F,| (16) 


Thus, Eq. (9) is reduced, within the approximations 
stated, to the nonhomogeneous linear partial differential 
equation 


or? dr \r dz? 
ra) wr fa) F, ra) - 
(« “) (aor) + f(r,z) (17) 


where nondimensional velocities u/wo, are de- 
noted i, %; when the blade forces’are prescribed, is, 
known from Eq. (13). 

The axial and tangential velocities corresponding to, 
this radial velocity follow from the equation of con- 
tinuity and the second equation of motion, respec- 
tively. Through integration of Eq. (8) the axial 
velocity is 

(w/w) -1=®@—1= (0/rdr)(ru)dB (18) 
while, by the appropriate approximation that lé — i 
<< 1, the tangential velocity is, from Eq. (11), 

— = —f u(0/rdr)(dr)dB 


The complete mathematical problem of calcylating 
the linearized radial, tangential, and axial velocity 
components in a turbomachine with prescribed blade. 
loading is given by Eqs. (13), (17), (18), and (19), to-. 
gether with the following boundary conditions: 


(19). 


s= 


+ | ®= 1; z= 
u=0; 
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Fic. 3. Jump of tangential velocity across the discontinuity 


and either 
F;(r, 2), F,(r, 2) or 


prescribed throughout the annular space occupied by 
the fluid. 


do(r, 2), F,(r, (21) 


LINEARIZED SOLUTIONS FOR THE VELOCITY 
COMPONENTS 


The linearized solution for the radial velocity com- 
ponent will be obtained by finding the appropriate 
Green’s function—that is, a solution G(r, 2; a, 8) which 
gives the radial velocity, consistent with the boundary 
conditions, induced at any point of a circle 7, z by a 
unit change in tangential vorticity of a vortex ring of 
radius @ at an axial coordinate 8. The complete radial 
velocity is simply the sum of such solutions correspond- 
ing to all changes of tangential vorticity. From the 
manner in which the problem was linearized, it is clear 
that the solution G(r, z; a, 8) may also be interpreted 
as the velocity induced at a point r, z by a cylindrical 
surface made up of annular vortices with radius a 
extending from the axial coordinate 8 to ~. It will 
be of particular interest to find the Green’s function by 
considering the case (Fig. 3) where all forces.are con- 
centrated in a plane z = 8, so that the mean tangential 


velocity changes discontinuously from to wu. 


This plane corresponds to a sheet of radial and tangen- 
tial vorticity normal to the turbomachine axis. Clearly, 
Of(r, z)/0z = 0 when z ¥ 8; Eq. (17) becomes homo- 
geneous and may be solved directly, using the given 
boundary conditions, for independent solutions on 
either side of the discontinuity. Now it is clear from 
the left side of Eq. (17) and the symmetry of the bound- 
ary conditions that any solution of Eq. (17) corre- 
sponding to a unit disturbance at z = 6 must be sym- 
metrical in z about 8. Therefore, both radial and axial 
velocity components are continuous at z = 6. If the 
notation [ f(r, z)], is used to denote lim,,» (f(r, 8 + €) — 
S(r, 8 — ©)), the jump in the value of f(r, z) across the 
discontinuity at z = 8; then the jump in tangential 
vorticity may be written as 


Bil = w (f(r, 2) 


1 > (22) 
= - or) (vor) + 


inasmuch as the axial velocity is continuous at z = £. 
Thus the continuity of the radial velocity, together with 


the prescribed discontinuity in tangential vorticity at 
z = B [Eq. (22)], is sufficient to join the independent 
solutions at z = B.- 

If the applied forces, and hence the discontinuity 
in tangential vorticity, be limited to the neighborhood 
of a circle ry = a and vanish for other radii, the solution 
obtained corresponds to the desired Green’s function 
or unit solution. If, on the other hand, there exists 
a finite jump of tangential vorticity for all radii, the 
solution bears the same relation to that for a continuous 
axial distribution of force as the lifting line solution for 
finite wings does to the solution for continuous chord- 
wise load distribution. 

If a and #@ are the nondimensional radial veloci- 
ties upstream and downstream of the discontinuity at 
z = 8, this new problem may be formulated as 


Ps) (*) 
or? or\ r oz? (23) 
with the boundary conditions 
a) = 0; z= 
=0; s=+0 (24) 
a” == 0; = 11, 72, t= 2 


and the auxiliary conditions at the discontinuity 


B) = B) 
ou ou 
oz oz L (95) 


A solution to Eq. (23) of the form 
> Uy(enr)(A,e* + B,,e~**) (26) 
n=1 


holds on each side of the discontinuity, where U;(e,7) is 
the linear combination of Bessel functions of order one 


= Jilenr) Valenti) — Silents) Vilenr) (27) 


and -the characteristic values e, are the roots of the 
transcendental equation 


Ui (€nt2) Yi (€n71) > Ji(€n?1) = O (28) 


Approximate values of the first six roots of Eq. (28) 
are tabulated in Jahnke-Emde.'® The conditions (25) 
relate the solutions #7 and a and, consequently, 
A,™, A,®, B,™, and B,®. 


= (29) 
Uilenr (By? — B,) = 1/2[ f(r, 
n=1 
In order to satisfy the conditions that the radial veloc- 


ity vanish at large distances upstream and downstream 
of the discontinuity (boundary condition 24 at + ~), 
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A, = B,® = 0. Hence, by applying the ortho- 
gonality relations of the Bessel functions*to Eq. (29), 
the coefficients are 


A, —e wv [ fla, da 
( * aU (Ena) 
B® = —¢ [ f(a, da 
ri 2v,,7€ En 


where a is introduced as the variable of integration and 
y, is the norm of the function U;(e,7) over the interval 


11, 


al *(€,a)da = 


2 


(31) 


The corresponding solution for the axial velocity fol- 
lows from Eq. (18) 


(32) 


OF A FLUID, THROUGH AN AX 


CIAL TURBOMACHINE 477 


where 


Uo(enr) = Solent) — Yo(enr) (33) 
and the functions y“(r) are determined by the initial 
distribution of axial velocity and the continuity of 
axial velocity at z = 8. Thus 


and applying these relations 
= 0 
yr) = Uolesr) X 
35 
JI, 


The linearized solution to the tangential velocity dis- 
tribution is found through integrating the known value 
of a” in Eq. (19). 


Substituting the known values for the constants and interchanging the order of summation, the complete 


solution for the single discontinuity becomes 


= 


where 


An examination of these solutions shows that they con- 


sist essentiallyin unit solutions of weight [f(a, 


summed over the radius from 7; to 72; the expressions in 
brackets represent solutions that give the velocities 
consistent with the boundary conditions induced at a 
circle r, 2 by a unit disturbance on the circle of radius 
a at the coordinate z = 8. These functions are de- 


noted 


> aU (enr) Ui(ena) 


n=1 


— 


G(r, 2; a, B) 


(39) 
aU o(€nr) Ui (ene) en 


H(r, 2; a, B) 


gas? Rape 


[2 > (1 


n=1 2€nVn 


= da (37) 


aU; Ui (ena) (38) 


da dz 


K(r,s; a, 8) = 2 (39 
cont.) 
(1 — 


The discontinuous approximations [Eqs. (36), (37), and 
(38)] are solutions of a linear problem and may be 
superposed with the convention that [ f(a, 2)]g at each 
value of 8 be interpreted not as the sum of the super- 
posed discontinuities but as the discontinuity corre- 
sponding to the sum of the individual tangential veloci- 
ties [see Eq. (22)]. This is necessary inasmuch as 
the approximate tangential velocity enters nonlinearly 
in evaluating the jump across the discontinuity. The 
solution for a continuous blade loading follows either 
by considering the limit of such a superposed sum of 
discontinuous jumps in which [f(a, 2)]g becomes 
[df(a, B)/2B]dB8 or by considering the original equa- 
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tions, Eqs. (17), (18), and (19), in relation to the unit 
solutions, Eqs. (39), which have been found. Hence, 
for any continuous loading, 


i = Ofla, B) G(r, z; «, B)\da dB (40) 
ri 
®-1= z; a, B)da dB + 


K(r, 2; a, B)\dadB (41) 
— = 


f G(r, 2; a, B)da dB dz 


(42) 


SPECIAL FORMS OF THE BLADE LOADING 


The case where the blade loading is similar at all 
radii and the effect of the radial force component may 
be neglected allows simplification of the expressions 
for the velocity components. Let R and Z be non- 
dimensional functions of r and z alone and jc be a 
parameter indicative of the total circulation about a 
blade of chord ¢. The radial and axial vorticities and 
the tangential velocity may be expressed in the form 


= —fcRZ’ 
= fc[R’ + (R/n)]Z (43) 
= 


where 


(€n7) 
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ent 


n=1 n 


n=1 


— = 2 Uilew) (Rr)’ T, )dz 


Wo 


T, 


/ 2 


where the prime denotes differentiation with respect to 
the appropriate variable. Then, according to Eq, 


(17), 


Using this relation in the solutions [Eqs. (40), (41), 
(42)] and interchanging the order of integration and 
summation, the following four fundamental integrals 


appear: 
= R(R'a + R)U(e,0)da 
(a@/r2)(R’a + R)Ui(e,a)da 


(45) 


The linearized approximation to the velocity compo- 
nents may then be separated into the velocity associated 
with a stationary blade row and the additional veloc- 
ity resulting from rotating a blade row of the same 
loading with an angular velocity w. Hence, 


ty + (wre2/ te 
—1=@, —1+4 (46) 
=) — + (wre/kc)ie 


) 


| 
| 
(47) 
| 
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Therefore, when the series 47 have been evaluated for a given type of blade loading, the velocities induced by 
any stationary or rotating blade row of similar loading 
that as the angular velocity is changed the blade loading 


follow directly from Eqs. (46). It should be noted 
and not the blade shape is preserved. 
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FLOW 


This flow may be approximated by a discontinuous 
solution or, more conveniently, by two discontinuous 
solutions—one associated with the velocities induced 
by the stationary blade row and the other with the ef- 
fect of rotation. If these two discontinuities are 
located at z = 8, and z = fz, respectively, the solutions 
follow from Eqs. (46) and (47) by the correspondence 


(2))2 — 
> 2 Ail 
Zo — Za 
0 
0 < Be 
7 (2) — 
Z 0 73 < 


where Z and Z® are defined by the approximate 
tangential velocity upstream and downstream of the 
blade row 

= 

= (49) 


and gc the circulation parameter remains constant as 
0. 
The values of 8; and @. are chosen to locate each of 
the discontinuities at the centroid of the loading it re- 
places. Thus, 6; and £; are solutions of the equations 


50) 
Sj, — 


eo. 


SINGLE Row oF STATIONARY OR ROTATING BLADES 


In the following example the velocities induced by a 
single stationary or rotating blade row are calculated, 
assuming that the radial force vanishes, using both a 
particular continuous chordwise loading and the corre- 
sponding discontinuous approximation. In addition 
to illustrating the general procedure of application, the 
numerical results provide a direct comparison of the 
continuous distribution and the discontinuous approxi- 
mation. 

Consider a single row of blade chord ¢ and prescribe 
the vorticity distribution &(r, z) as shown in Fig. 4. 


OF A FLUID THROUGH AN AXIAL TURBOMACHINE 
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Fic. 4. Distribution of load on the blade row. 


This distribution is similar at all radii and increases 
linearly from the blade root to tip. If the correspond- 
ing physical blade row is supposed to have m blades, 
the linear vorticity distribution of the blade at any 
radius 7 is given by 


E(r, 3) = (2Qer/n)k(r, (51) 
and the corresponding total circulation is 
T(r) = (34/2n) &c?(r/r1) (52) 
The tangential velocity becomes, from Eq. (13), 
% = 0 
[245] 
3 1 c 
and the parameters R and Z are 
R=r/n (54) 
Z = (e/e) + (1/2), $0 | 


Z= 3/4 


2 6/2 


The basic integrals, Eqs. (45), are then evaluated by 
straightforward calculation, and the results are given 
in the Appendix. 
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Fics. 5. Velocity component inducted by a blade of finite 
chord. Fic. 5a. Axial velocity distribution; stationary blade 
row. 


Fic. 5b. . Axial velocity distribution associated with angular 
velocity of blade row. 


Using these relations the velocity perturbations as- 
sociated with the stationary or rotating blade row fol- 
low directly from Eqs. (47). For the case where 
= 0.6, (r2 — m1) /¢ = 2, and the vorticity distribu- 
tion is of the form given in Fig. 4, the induced velocity 
distribution has been calculated for various distances 
from the center of the blade row. The results are pre- 
sented in Fig. 5. 

Figs. 5a and 5b indicate that, for the blade loading 
under consideration, an appreciable change in the axial 
velocity profile takes place both before the fluid enters 
and after it leaves the blade row, the coordinates 
z/¢ = +0.5 marking the termini of the blade row. 
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Fic. 5c. Tangential velocity distribution; stationary blade 


row. 
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Fic. 5d. Tangential velocity distribution associated with 


angular velocity of blade row. 
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FLOW OF A FLUID THROUGH AN AXIAL TURBOMACHINE 


Fic. 5e. Radial velocity distribution; stationary blade row. 


These axial velocity changes and the corresponding 
ones for the tangential velocity (Figs. 5c and 5d) are 
those involved in considering the interference of ad- 
jacent blade rows, and for the sake of simplicity it is 
convenient if they can be estimated from the discon- 
tinuous approximation. The radial velocity disturb- 
ances (Figs. 5e and 5f), which are of practical import- 
ance only since they produce the changes in tangential 
and axial velocity, exhibit a nonsymmetry about the 
centerline (z/c = 0) of the blade row which results from 
the unsymmetrical blade loading. The disturbance as- 
sociated with the stationary blade attains its maximum 
value behind the centerline, while the additional dis- 
turbance associated with the rotation attains its maxi- 
mum value ahead of the centerline. This difference is a 
direct result of the load distribution producing the flow 
in each case and is clearly reflected in the position of 
the discontinuities appropriate for approximating each. 
Substituting the known load distribution into Eqs. 
(50), it follows that 

Bi = 0.026c 
Bo 


which agrees with the observed nonsymmetry of the 
flow. 

The discontinuous approximation to the axial veloc- 
ity variation associated with the stationary blade row 
follows directly from the correspondence given by Eqs. 
(48) where the discontinuity is located at 6,/c = 0.026. 
The.axial velocity variation given by the discontinuous 
approximation is compared with the continuous solu- 
tion in Fig. 6 and clearly provides a reasonable approxi- 
mation except, of course, in the immediate vicinity of 
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Fic. 5f. Radial velocity distribution associated with angular 
velocity of blade row. 
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Fic. 6. Comparison of axial velocity distribution given by the 

discontinuous solution with that for finite chord. 
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the lifting line. This discrepancy near the lifting line 
is to be expected because of the excessively high radial 
velocity induced in this region which causes a corre- 
spondingly high axial gradient of the axial velocity per- 
turbation. 
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APPROXIMATION TO THE VELOCITY COMPONENTS 


The distribution of velocity far downstream of the blade rows, which is analogous to the downwash in the 
Trefftz plane for the theory of finite wings, may be calculated with comparative ease. For upon integration with 
respect to z, Eqs. (10) and (17), give, far downstream, 


(0/0r)(@.. = fir, ) — fir, (57) 


where @_., and @,, correspond to the distributions of axial velocity far upstream and downstream, respectively, 
of the blade row. Thus, 


— Wo = S (f(r, — fir, ))dr + 
so that if the approximate tangential velocity distribution is given by % = fcg(r, 2), it follows that 


(fc /wo)? /wo)? r or 


r20r 


(foc /w)? (fc 


(58) 


where the constants C, and C; are evaluated to make 
(@ — 1).r dr = 0, (@.).r dr = 0 (59) 


Now for the case in which the velocity far upstream is uniform, it follows from the general discontinuous ap- 
proximation to the axial velocity [Eqs. (46), (47), (48)] that 


lim {25 > 


n=1 


1 


rez — n(a \ = hy(r) + ri h2(r) (60) 


wet 1 2€nVn? 


This suggests the approximation 


and (61) 
13 < Be 


where § is the blade aspect ratio (7.2 — 7,)/c and the \ are undetermined constants. Similarly, using the continu- 
ity relation Eq. (8) and from Eq. (19), the corresponding radial and tangential velocity components are 


tty / (oc /wo)? = — (v2 — Ri | 
dels — Be (62) 
and 
r 
(63) 
(foc /wo)? Wo (r) r (e 1); Bo 
where 
rki, V/s rh, 2(r)dr +C (64) 


and the constant of integration is evaluated to satisfy the boundary conditions on the radial velocity. 
If the radial velocities [Eq. (62)] satisfied the homogeneous differential equation, Eq. (23), the velocity com- 
ponents Eqs. (61)—(63) would be true solutions of the linearized problem. Hence, the values of \; and i, will be 
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chosen so that the radial velocity satisfies Eq. (23) except for a minimum weighted mean square error over the 
entire field of flow. From the first of Eqs. (62) and Eq. (23), 


and the value of \; is found by minimizing the integral 


khi(r) [(s — 61) 


2 ’ 2 2 
r Fe 


By integrating with respect to 2, minimizing with 
respect to Ai, and simplifying the result through inte- 
gration by parts, the parameter ), is 


2d 
ky r dr 


where each of these integrals is positive and, hence, the 
values of the parameter are rea]. Similarly, for the 
solutions of the additional velocity associated with a 
rotating blade row, 


re 2 2 1/2 
r dr 


In the example for which detailed calculations were 


(66) 


== 


(67) 


carried out in the previous section, v°? = 0; ? = 
(3/4) fc(r/r1). It follows directly that 
al +1 9 
32L\n 
(68) 
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where 7;/r2 = 0.6. The appropriate values of the ex- 
ponents and are then 


and 


The exponential approximation to the velocity dis- 
tribution may be found very simply by applying Eqs. 
(61)-(63). In particular, the axial velocity perturba- 
tion associated with the stationary blade row is 


®, — 1 


(c/w)? ~ 


9| /r? r \? 
1 2 1.29(2 — c) 
ai («) |. 
(71) 
9 r \2 


and is compared in Fig. 7 with the general linearized 
solution for finite blade chord. This exponential ap- 
proximation appears, at least in this case, to be suf- 
ficient for estimating the general order of the velocity 
components. Their simplicity makes these relations 
useful in estimating the interference of adjacent blade 
rows. 

A further simplification, useful in the neighborhood 
of the blade row, follows by. approximating the exponen- 
tial function in Eqs. (61)—(63) in the form 


+A[(z—8)/c 


The axial velocity upstream of the blade row is given, 
for example, by 


,2< 
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Fic. 7. Comparison of axial velocity distribution given by ex- 
ponential approximation with that for finite chord. 
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and similar relations may easily be obtained for the other velocity components. 


Appendix 


VALUES OF THE BASIC INTEGRALS 


= — (r2/r1)?Uo(enr2) 
T, = (n/re)T,© 


The functions J’, must be evaluated in four separate cases corresponding to the four different representations 
of Z. 
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Theoretical Pressure Distributions for a Thin 
Airfoil Oscillating in Incompressible Flow 


E. E. POSTEL* ano E. L. LEPPERT, JR.* 
Lockheed Aircraft Corporation 


ABSTRACT 


Formulas are derived for the theoretical pressures acting on a 
thin airfoil and aileron combination oscillating with harmonic 
motion of small amplitude in incompressible flow. 

Calculations of pressure distribution are made at several 
representative values of the frequency parameter 1/% for airfoil 
translation, for airfoil rotation about the quarter-chord, and for 
hinged-surface flapping motions. The tabulated data thus ob- 
tained are plotted. , 


INTRODUCTION 


FORCES AND MOMENTS for the har- 
monically oscillating airfoil have previously been 
derived by Theodorsen' without the pressure formulas 
having been explicitly stated. A number of German 
papers?~® have been written covering various phases of 
the subject. Of these, Dietze? has developed pressure 
formulas (for translation, pitch about the leading edge, 
and control surface flapping) which, when referred to 
the quarter-chord axis, can be shown to be in agreement 
with those developed in this paper. In these German 
investigations, as in Theodorsen’s, the primary objec- 
tive was to obtain the lifts and moments; therefore, 
computations were not made to show typical pressure 
distributions. 

In a recent investigation to determine the magnitude 
and phase of aileron hinge moments due to an aileron 
oscillating harmonically in a transonic airstream,’ 
pressure-measuring units were installed in the wing 
and aileron upper and lower surfaces to measure the 
oscillating pressures. The pressure pickups were de- 
signed specifically for the conditions existent in that 
test and were, therefore, not satisfactory for measuring 
the oscillating pressures at the lower dynamic pressures 
of normal flight speeds. Formulas were derived and 
computations for theoretical oscillating ptessures were 
made for the range of parameters expected in that test 
in order that comparisons with the experimental pres- 
sures could be made. The derivations and the tables 
of computed pressures for the incompressible flow are 
presented below in anticipation of the development of 
pressure-measuring techniques to check the theoretical 
flutter coefficients in the usual velocity range. 


NOMENCLATURE 


The nomenclature used in the following developments is that of 
Theodorsen,! with the exception of definition of several additional 
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quantities useful in the formation of the pressure formulas: 


h = vertical translation of airfoil 

a@ = airfoil pitch (rad.) 

B = aileron rotation (rad.) 

# = distance aft of mid-chord measured in _half-chord 
lengths 

Xo = x position of vortex elements measured in half-chord 
lengths 

a = distance from mid-chord to axis of rotation in half- 
chord lengths 

c = distance from mid-chord to axis of flapping surface 
hinge line in half-chord lengths 

b = half chord 

¢e = velocity potential due to sources and sinks (noncircu- 


latory flow). Note: Subscripts ahafB define the 
¢. for the motion described by the subscripts 


¥r = velocity potential due to vorticity (circulatory flow) 

p = total pressure difference between the upper and lower 
surface 

De = pressure difference due to noncirculatory flow 

pr = pressure difference due to circulatory flow 

p = density 

v = free-stream velocity 

w = frequency of oscillation 

k = bw/v reduced frequency parameter 

= vorticity 

= density of vorticity AT 


_ U dx. 
= Udxo 
—1 
Xo 
Cw = F+iG 


+1 
Vx? — 1 
i = -1 
N 
To = V1 — c + ¢ 
(1 — 2c) + W1 — (2 0) 


e dxy 


Tu 
The additional terms defined are: 
Qe = — V1 — (d¢/dx) 
Qa) = Qa atx = 1 
p 


DERIVATION OF PRESSURE FORMULA IN TERMS OF 


VELOCITY POTENTIALS 


The mathematical development used by Theodorsen 
in obtaining the forces and moments on an oscillating 
surface defines the velocity potential for noncirculatory 
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flows and for those due to the surface of discéntinuity 
behind the wing. In the following, the “noncircula- 
tory’’ pressure is stated in terms of the noncirculatory 
potential, and the circulatory pressure is stated in terms 
of the circulatory potential according to Theodorsen. 
The total pressure difference is then stated in terms of 
functions of the noncirculatory potentials by connecting 
the circulatory and noncirculatory potentials through 
the Kutta condition. 


p = bh + br 

pe = 4 


1 xo? — 1 —1 
Ud 
| 


x02 — 1 (x0°— x) 


2rb 
| 1 
V1 — x? 
xX +x 


— 2pv 1 ) 
(7; — x Vi? — 1 @) 


The magnitude of the circulation is determined by 
the Kutta condition, which requires that no infinite 
velocities exist at the trailing edge, or, at x = 1, 


© (xp? — 1) + 1 — x? 
Vx? — 1 (xo — x) U 


the total pressure difference is thus 
Qa) 1 0¢, 


={- + + CQa (7) 


CALCULATION OF Q’S IN TERMS OF ¢& 


The noncirculatory velocity potentials for the various 
motions of the airfoil-aileron combination may all be 
expressed in terms of four basic functions. Accord- 
ingly, values of Qa, and + xQua)) are ob- 
tained below: 

V1 — x? 
Qe = (—V1 — x4)(OV1 — x*/da) = x 


= 1 
+ = 0 


For = — x? 


For = 


fa) 
(er + + a t+ + + = finite 


Introducing the definition of 
Qe) = —V1 — x? (d¢/dx) 


at x = 1, and for gp from O¢gy,/Ox according to Theo- 
dorsen’s Eq. (VII) 


Qay = (1/22) x0 1/V x0 — 1) Udxo (3) 


or 
Qa (1/2r) fi [(xo + 1)/V x0? 1] U dxo (3a) 
then 


1 V x9? 
1 ws 
on 1 U Xo 
= — /V (xo? — 1) 1) U dxo (4) 


subtracting Eq. (4) from Eq. (3a) 
Qw = CQ a) = (1/2r) f° U dxo/V x02 — (5) 
stating Eq. (2) in terms of Eqs. (4) and (5) 


= + CQa(1 x) ] (6) 


= + CQm(1 — 


Qe) = —(1 — 2x*) 
Qa = 1 
+ = (1 — x)(1 + 2x) 


For ¢3 = (x — c) log N 
= V1 — x? log N 
=-Vi-@ 
—Qe + = (1 x)V1—2+ V1 —x?log N 


For g = (x — c) log N 
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Qw = — — — c) V1 — log N | — Qu + xQu = + 
= 2 (x — c)V1 x? log N 


PRESSURE EQUATIONS 


The pressure equations are now obtained by writing the total velocity potential for the motion desired and sub- 
stituting values from the four basic ¢ functions in Eq. (7). 
Thus, for # motion, assuming simple harmonic motion of translation, 


ye = hbV1 — x? = iwhbV1 — x? 
10 1 
=-— + + + CQa (1 — x)] 
— x? 


+ 2+ Cy (1 x)] 


= (OV 1 — x?) (tw) + 


—pb b 1 
1—x 


2pv7he 


— x? + ikCV (1 — x)/(1 + x) 


or 


—p/2pu*bhe = —V1 — x? + (1 — +x) 


Cy = —V1 = x? + — = + (9) 


Similarly, but in an algebraically somewhat more complicated manner, the pressure equations for a and 8 mo- 
tions are obtained: 


a Motion (about quarter-chord, a = —1/2) 
—) = 1 1 1-<x —-x 
(: - 5 Vi + l+x +2 (10) 
separated into its real and imagina1 y components 

C, (real — sai 
> (real) = Vi (11) 

C, (imag.) = — x2+ + (05 +x + (12) 
k k l+<x 


Note that in the above equations the rotation is taken about the quarter-chord axis. This leads to a simple 
variation for several of the force coefficients. Rotation about any desired axis can be obtained by combining a 
and / motions. 


8 Motion (rotation about hinge line at c) 
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1 (x — 2c) 1 ~ 


2a 
[z= Tut | fi 
— — 14) 
wk? mk INJ + % 


cos~!¢ 2 
G (imag.) = Vi — — 0) log N+ 


wk 
Es TuF | cos~'c¢ (1 + 2x) (15) 


The pressure equations are readily integrated into the real and imaginary coefficients of Theodorsen. The rela- 
tion between them is stated below: 
(2/r) £-t" C, dx = R + Ii (force) 
(2/m) + a) dx = R+ Ji (total moment about axis at a) 
(2/n) C(x dx =R+Ki (flapping surface moment about hinge line) 


TABLE 1 
Pressure Coefficients for 4 Motion 


Station 
Per Cent ———1/k = ——1/k = 2.94—— ——l/k = 3.33——~ ——l/k = 4.17—— 
Chord Real Imag. Real Imag. Real Imag. Real Imag. 
0 0 0 . 
5 —0.436 1.792 8.293 2.169 9.662 2.946 12.694 
10 —0.600 0.934 5.708 1.193 6.650 1.727 8.736 
15 —0.714 0.503 4.529 0.709 5.277 1.133 6.932 
20 —0.800 0.222 3.805 0.395 4.413 0.752 5.824 
25 — 0.866 0.089 3.295 0:169 3.839 0.478 5.044 
30 —0.917 —0.136 2.906 —0.004 3.386 0.269 4.448 
35 — 0.954 —0.257 2.593 —0.139 3.021 0.103 3.969 
40 — 0.980 —0.354 2.330 —0.248 2.715 —0.030 3.567 
45 —0.995 —0.430 2.103 —0.334 2.451 —0.137 3.219 
50 — 1.000 —0.489 1.903 —0.402 2.217 —0.224 2.912 
55 —0.995 —0.533 1.721 —0.454 2.005 —0.298 2.634 
60 —0.980 —0.562 1.554 —0.492 1.810 —0.346 2.378 
65 —0.954 —0.579 1.396 —0.515 1.627 —0.384 2.137 
70 —0.917 —0.582 1.246 —0.525 1.451 —0.409 1.906 
75 —0.866- —0.571 1.099 —0.521 1.280 —0.418 1.681 
80 —0.800 —0.544 0.951 —0.501 1.108 —0.412 1.456 
85 —0.714 —0.499 0.799 — 0.463 0.931 —0.388 1.223 
90 —(.600 —0.430 0.634 —0.401 0.739 —0.341 0.971 
95 —(0. 436 —0.319 0.437 —0.298 0.509 — 0.258 0.668 
100 0 0 0 0 0 0 0 
RESULTS 
Pressure coefficients, calculated in accordance with 
the derived formulas of Eqs. (9), (11), (12), (14), ana 
(15), are computed for airfoil translation, rotation 
about the quarter-chord, and control surface flapping. 
— c + in Tables 1, 2, and 3, respectively. The computations 
= 1 shiek ies ' are made at 1/k values of 4.17, 3.33, 2.94, and 0, for 5. 
” | 4 } = per cent chord intervals for h and a motions. Compu- 
h tations for 8 motion are made at 1/k values of 4.17, 
~< 3.33, and 2.94 for hinge line at 75 per cent chord and 
at 1/k = 2.94 for hinge line at 50 per cent chord. 
The real and imaginary components are plotted 
separately in Figs. 2 through 7. The real component 
is the value of the pressure coefficient at the time the 


Fic. 1. Airfoil parameters. displacement is maximum. The imaginary component. 
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TABLE 2 
Pressure Coefficients for a Motion, a = —0.5 
Station 
Per Cent ——l1/k =0 ~ 1/k = 2.94——~ ——1/k = 3.33—~ ——1/k = 4.17——— 
Chord Real Imag. Real Imag Real Imag Real Imag 
5 0.022 26.598 —2.106 34.793 —3.381 56.247 —6.845 
10 0.060 18.261 0.315 23.901 —0.327 38.667 —2.211 
15 0.107 14.430 1.650 18.906 1.328 30.622 0.229 
20 0.160 12.054 2.563 15.814 2.449 25.658 1.859 
25 0.217 10.361 3.238 13.617 3.275 22.143 3.053 
30 0.275 9.054 3.755 11.925 3.907 19.444 3.966 
35 0.334 7.989 4.151 10.551 4.394 17.258 4.674 
40 0.392 7.088 4.451 9.390 4.766 15.418 5.221 
45 0.448 6.304 4.668 8.382 5.039 13.824 5.634 
50 0.500 5.607 4.811 7.487 5.224 12.409 5.930 
55 0.547 4.977 4.884 6.677 5.329 11.129 6.117 
60 0.588 -4.398 4.889 5.934 5.354 9.952 6.202 
65 0.620 3.861 4.825 5.241 5.301 8.853 6.186 
70 0.642 3.356 4.690 4.587 5.166 7.809 6.064 
75 0.652 2.876 4.476 3.962 4.941 6.803 5.829 
80 0.640 2.414 4.170 3.354 4.612 5.815 5.465 
85 0.607 1.958 3.751 2.748 4.155 4.816 4.941 
90 0.540 1.496 3.172 2.122 3.519 3.763 4.199 
95 0.414 0.987 2.318 1.418 2.575 2.547 3.081 
100 0 0 0 0 0 0 0 
TABLE 3 
Pressure Coefficients for 8 Motion 
Station c=0 c= +0.5 
Per Cent ——l1/k = 2.94—~ _ ——1/k = 1.0——~ 1/k = 2.94 ——1/k = 3.33—~ ——1/k = 4.17— 
Chord Real Imag Real Imag. Real Imag Real Imag Real Imag. 
21/2 7.558 —4.949 10.397 —6.627 18.283 — 10.924 
5 10.472 —3.043 0.400 —0.723 5.650 —3.373 7.739 —4.506 13.516 —7.511 
10 8.205 —2.149 0.337 —0.457 4.453 —2.198 6.052 —2.985 10.439 —4,999 
15 7.440 —1.705 0.323 —0.326 4.041 — 1.639 5.455 —2.246 9.305 —3.814 
20 7.182 —1.387 0.325 —0.247 3.883 —1.297 5.212 —1.796 8.802 —3.091 
25 7.205 —1.110 0.335 —0.172 3.851 —1.002 5.141 —1.419 8.605 —2.505 
30 7.447 —0.841 0.350 —0.116 3.895 —0.777 5.175 —1.130 8.592 — 2.058 
35 7.930 —0.555 0.371 —0.067 4.012 —0.588 5.308 —0.888 8.748 — 1.687 
40 8.776 —(.232 0.397 —0.021 4.182 —0.412 5.512 —0.667 9.023 —1.351 
45 10.487 0.168 0.429 0.024 4.418 —0.246 5.803 —0.459 9.440 —1.042 
50 ey 0.805 0.469 0.069 4.731 —0.082 6.194 —0.256 10.017 —0.742 
55 10.292 1.422 0.521 0.118 5.154 0.087 6.727 —0.049 10.821 —0.443 
60 7.815 1.756 0.591 0.171 5.744 0.270 7.476 0.173 11.964 —0.129 
65 6.457 1.967 0.697 0.235 6.645 0.481 8.625 0.425 13.735 0.221 
70 5.418 2.085 0.904 0.326 8.446 0.772 10.931 0.768 17.316 0.688 
75 4.546 2.122 © 0.475 © 1.231 © 1.236 © 1.381 
80 3.770 2.078 0.813 0.612 7.748 1.651 10.025 1.789 15.861 2.025 
85 3.044 1.944 0.544 0.636 5.481 1.751 7.107 1.916 11.283 © 2.127 
90 2.326 1.698 0.362 0.592 3.921 1.646 5.096 1.811 8.117 2.123 
95 1.543 1.275 0.210 0.467 2.501 1.308 3.259 1.444 5.208 1.7 
100 0 0 0 0 0 0 0 0 0 0 


is the value as the displacement passes through zero. 
When the imaginary component is positive, the pressure 
exerts a damping force. 


In order to show more readily the relative phase and 


magnitude of the pressures for variation in chord posi- 
tion, pressure coefficients at three values of 1/k are 
plotted in complex form for h motion in Fig. 8. Fig. 9 
depicts the change in phase and magnitude of 8 motion 


CONCLUSIONS 


pressures due to change in hinge line position. The 
value of 1/k is 2.94. 


These formula and curves make possible a comparis 
with experimentally determined pressures which cal 
point the way to empirical modification of the two 
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Fic. 8. Complex pressure coefficient, # motion. 


dimensional coefficients and may make possible a 
modification of theory. Of equal importance is the 
possibility of comparing three-dimensional experi- 
mentally determined pressures with present strip- 
theory methods of estimating three-dimensional effects. 
Development of instrumentation has progressed suf- 
ficiently to indicate that pressure measurements are 


* Dietze, The Air Forces on a Harmonically Oscillating Self. 
Deforming Plate, Luftfahrtforschung, Vol. 16, No. 2, p. 84, 1939, 

3 Jaeckel, On the Calculation of the Circulation Distributions fora 
Two-Dimensional Wing Having Periodic Motion, Luftfahrt- 
forschung, Vol. 16, No. 3, p. 135, 1939. 

4 Schwarz, Calculation of the Pressure Distribution of a Wing 
Harmonically Oscillating in Two-Dimensional Flow, Luftfahrt- 
forschung, Vol. 17, Nos. 11, 12, p. 379, 1940. 

5 Sohngen, Determination of the Lift Distribution for Optionally 
Non-Uniform Movement, Luftfahrtforschung, Vol. 17, Nos. 11 
12, p. 401, 1940. 

6 Dietze, Law of Aerodynamic Force of a Jointed Plate in Har- 
monic Motion, Luftfahrtforschung, Vol. 18, No. 4, p. 135, 1941. 

7 Lockheed Report No. 5808 (unpublished). 
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Stability of Boundary Layers and of Flow 


Entrance Section of a Channel | 


ELIZABETH HAHNEMAN,* J. C. FREEMAN,? M. FINSTON? 


Brown Unwersity 


SUMMARY 


Lin’s approximate formula for the minimum critical Reynolds 
Number is used to study the entrance flow into a channel, the 
boundary layer with pressure gradient, and the boundary layer 
with suction. It is found that the results thus obtained compare 


favorably with much more elaborate calculations wherever they 


areavailable. In the case of the entrance flow into a channel, the 
new result is found—that the entrance section is comparatively 
more stable than the fully developed parabolic velocity distribu- 
tion. 


INTRODUCTION 


T NHE STABILITY CHARACTERISTICS of several standard 


types of flows have been investigated by various 
authors. Calculations were also carried out in certain 
cases to determine the effect of various factors, such as 
pressure gradient and suction upon the stability of the 
boundary layer. Such results are no doubt of interest, 
but in order to obtain them a series of elaborate calcu- 
lations has to be carried out. 

In studying the general problem of hydrodynamic 
stability of an incompressible fluid, Lin' gave some 
simple rules for the calculation of the critical Reynolds 
Number when the velocity distribution is known. These 
rules were later generalized by Lees? to the case of the 
boundary layer in a compressible fluid. They are par- 
ticularly suitable for the problem at hand, since they 
readily give the critical Reynolds Number once the 
steady flow is known. 

In the present paper, several cases are considered by 
the approximate formulas of Lin. The results are 
compared with existing calculations whenever they are 
available. In all such cases, the comparison is ex- 
tremely satisfactory. Some are presumably better 
than existing ones obtained from more elaborate calcu- 
lations. 

The cases considered are: (1) the stability of flow 
through a two-dimensional channel in its entrance sec- 
tion; (2) the influence of pressure gradient on the sta- 
bility of the boundary layer; and (3) the influence of 
suction on the stability of the boundary layer. 

In the first case the critical Reynolds Number is 
given by 

R = (30u;'/c*) [(wi'/c) w'dy]” (1) 
where c is the value of w(y), with y determined by 
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—rw,'(1 — 2y)(ww"/w’*) = 0.58 (2) 
and 
y = (wi'y/w) - 1 (3) 


In these formulas, w(y) is the velocity distribution; 
y1 is the coordinate of the solid boundary; ye is the 
coordinate for maximum velocity, dashes denote dif- 
ferentiation with respect to y, w,’ = w’(y), ete.; and 
R is the critical Reynolds Number. Note that all 
these quantities are dimensionless and the formulas 
hold for any consistent system of reference quantities 
chosen. 

For the boundary-layer calculations [cases (2) and 
(3) ] the formula for the critical Reynolds Number is 


R = 25w,'/c* (4) 


where the quantities w,’ and c have meanings as above 
but the Reynolds Number must be defined with respect 
to the free-stream velocity, although any reference 
length may be chosen. Usually, the displacement 
thickness is used. 

These calculations were made separately at the sug- 
gestion of Prof. C. C. Lin. The authors are indebted 
to him for close supervision of the several calculations 
and coordination of the results. 


STABILITY OF FLow THROUGH A CHANNEL 


The velocity distribution in the entrance section of 
the channel changes gradually from a uniform one at 
the entrance into a parabolic one at large distances from 


R 
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v 


Fic. 1. Curve showing the dependence of the minimum 
re Number on the distance from the entrance of a 
channel. 
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the entrance. The whole series of distributions was 
calculated by Schlichting. A detailed discussion of 
his results is given elsewhere.* These calculations were 
accurate for the velocity, although the derivatives of 
velocity cannot bé obtained with great accuracy. How- 
ever, despite this inaccuracy, a consistent curve can be 
obtained for the dependence of the critical Reynolds 
Number on the inlet distance. Fig. 1 shows such re- 
sults. The ordinate is the Reynolds Number defined 
in terms of the maximum velocity for the parabolic 
section and the half width of the channel. It is seen 
that the value decreases monotonically from infinity at 
the entrance to the value for parabolic distribution at 
large distances. It is clear that the entrance part may 
be treated as a boundary layer and, hence, has a finite 
critical Reynolds Number based on its thickness. Since 
the Reynolds Number plotted is based on the channel 
width, it becomes infinite. 

This result is significant in that any instability in a 
carefully controlled experiment would occur in the re- 
gion of parabolic velocity distribution, thus making the 
experimental check possible for the stability of such 
velocity distributions. The entrance section should 
not introduce serious complications. 


EFFECT OF PRESSURE GRADIENT ON THE STABILITY OF 
THE BOUNDARY LAYER 


As has been indicated before, the calculation of sta- 
bility depends on the velocity distribution w(y). Ap- 
proximate forms for w(y) may be found by various 
methods. Asis now well known, w” plays an important 
role in the stability problem. Schlichting’ modified 
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Fic. 2. Values of R, for Hartree’s profiles and Schlichting’s P, 
; profiles obtained from Lin’s formula. 
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Fic. 3. R. for Schlichting Ps profiles. 


the Pohlhausen method for the determination of w, so 
that a more accurate calculation of the second deriva- 
tive may be made, by taking a sixth order polynomial 
and introducing the w’’’ into the boundary conditions 
to determine the coefficients. For the Blasius profile, 
where the exact solution is known, Schlichting’s solu- 
tion shows a marked improvement over the Pohlhausen 
fourth order polynomial solution. 

The form of w is then given by 


w= F(y) + dG(y) (5) 
where 
F(y) = 2y — + — 
G(y) = — + yt — + 


and \ is the well-known Pohlhausen parameter de- 
pending on the pressure distribution.* Schlichting 
has tabulated F, G, F’, G’, F”, and G” for y, the bound- 
ary-layer coordinate, equal to 0, 0.05, 0.10, ... 1.000. 
These values and Eqs. (2), (3), and (4) are used to cal- 
culate the results shown in Fig. 2, noting that the re- 
sulting value of R is based on the boundary-layer thick- 
ness and must be multiplied by the factor 


6/86 = (2/7) — (A/105) 
to yield R based on the displacement thickness. RK is 
plotted against if in Fig. 2. Separation occurs at 


* In terms of the pressure gradient 0p/0x, the boundary layer 
thickness 6, the free-stream velocity w., the density of the fluid 
p, and the kinematic viscosity v, the parameter \ is defined by 


A = — (6*/pw*v)/(0p/ 0x) 
+ \, is a parameter based on 6,, the displacement thickness. 


AL = (6,?/vpw?) (Op/ Ox) 


and for the Schlichting profiles 
Ai = A[(2/7) — (A/105)]}? 
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TABLE 1 
Minimum Critical Reynolds Numbers for Velocity Profiles Re- 
sulting When the Suction Velocity », ~ 1/*/x on a Flat Plate 


R 
c (Ulrich) (Lin’s Formula) 
0 (Blasius Case) 575 465 
0.5 2,986 2,480 
1.0 9,550 8,510 
1.5 19,100 13,200 
TABLE 2 


Minimum Critical Reynolds Numbers for Velocity Profiles Re- 
sulting When the Suction Velocity v, = Constant on a Flat Plate 


& R (Ulrich) (Lin’s Formula) 
0.005 1,122 870 
0.02 1,820 1,460 
0.08 3,935 3,700 
0.18 7,590 6,750 
0.32 13,500 10,600 
0.50 21,900 13,600 
oo 70,000 (Bussman & Munz) 40,000 
50,000 (Pretsch) 
*)} = —10, for which the approximate rule gives 
R=0. 


The results obtained are in good agreement with 
Schlichting’s more elaborate calculations (see Fig. 3). 
For the Blasius profile, Schlichting gave the value R = 
575 for the exact velocity distribution and 645 for the 
polynomial approximation. Lin’s revised calculation 
for the exact velocity distribution gives the value 420. 
The approximate rule gives the value 502 for the exact 
velocity distribution and the value 558 for Schlichting’s 
polynomial approximation. 

The stability of profiles considered by Hartree was 
also examined. It was found (Fig. 2) that the critical 
Reynold’s Number for a given \ was within 20 per 
cent of that found by using the Schlichting-Ulrich ap- 
proximation [Eq. (5)]. Since the Hartree profiles are 
exact, the extent of agreement between the two calcu- 
lations indicates that the Schlichting-Ulrich approxima- 
tion gives good values for w”. It should be noted that 
only first differences were used in calculating the results 
based on Hartree’s work.*® 


STABILITY OF VELOCITY PROFILES RESULTING FROM 
SucTIGN ALONG A FLAT PLATE 


The approximate formula developed by Lin! is now 
applied to calculate the minimum critical Reynolds 
Number of the velocity profiles resulting from suction 
along a flat plate. The stability of some of these pro- 
files has been calculated by Ulrich,® 7 and the results 
are compared here. 

Ulrich® calculated the stability of the profiles result- 
ing from a suction velocity 1, ~ 1/+/x at the plate 
(where x is measured along the plate) by the elaborate 
methods of Tollmien and Schlichting. The profiles 
have been computed by Schlichting and Bussman.* 
They vary with the parameter. 


C = co(wel/v)'”* (6) 


where the suction coefficient 
= 2Q/(pw,*bl) 


and w» is the free-stream Velocity; /, the length of plate 
considered; v, the viscosity coefficient; Q, the total 
mass of air sucked in per unit time; p, the density of 
the air; and ), the width of plate considered. Ulrich’s 
values of the minimum critical Reynolds Number for 
various values of C are given in Table 1. 


The calculation of the critical Reynolds Number by 
Lin’s formula was performed on the original profiles 
from the accurate calculations of Schlichting and Buss- 
man. Thus one source of error was avoided. These 
results are also included in Table 1 and may be com- 
pared with those of Ulrich. 


Ulrich’ also calculated the stability of the profiles at 
various distances along a plate with a constant velocity 
of suction. According to Ulrich the profiles for this 
case were computed by Iglisch'! and vary with the 
parameter £ given by 


é = 0,71 /wev (7) 


The minimum critical Reynolds Number for various 
values of £, as computed by Ulrich, are shown in Table 
9 

Lin’s! formula, Eqs. (2), (3), and (4), was used on the 
Iglisch profiles, and these results are also in Table 2. 
The values obtained by Lin’s formula are 15 to 30 per 
cent lower than those obtained by the more laborious 
calculation of Ulrich, which is good agreement for sta- 
bility calculations. 

All of Ulrich’s calculations involve a curve given by 
Tollmien and Schlichting which is incorrect. Correct 
values of ,this curve, F(z), have been given by Lin,! 
and a nearly correct curve is used by Pretsch® in making 
his stability calculations. It is interesting to see how 
the use of a correct curve affects the value of the mini- 
mum critical Reynolds Number. For the Blasius case, 
using the incorrect curve, Schlichting obtained R = 
575; using the correct curve, Lin got R = 420. For 
the case of asymptotic suction, using the incorrect curve, 
Bussman and Miinz" obtained R = 70,000, and Pretsch, 
using a nearly correct curve, got R = 55,000. 

Thus it can be surmised that, since the actual values 
of R in the cases £ = 0 and £ = © are lower than those 
obtained by the methods used by Ulrich, Lin’s formula 
gave even more accurate results than indicated in 
Tables 1 and 2. 
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Letter + 


Dear Sir: 

In Professor Goldstein’s Wright Brothers Memorial Lecture 
(“Low Drag and Suction Airfoils,”’ JouRNAL OF THE AERONAUTI- 
CAL Science, April, 1948) a hint was made (page 209, lower 
right-hand column) that more than a direct decrease in drag fol- 
lowed from the drawing off of the boundary layer and re-energiz- 
ing it. Perhaps a more elementary demonstration of this phe- 
nomenon may be shown as follows: 

Assume a wing in a stream at velocity V, and further imagine 
for the purpose of simplicity that the boundary layer is not com- 
posed of a quantity of air Q, which is slowed to RV where k < 
1.0, but of a smaller quantity Q,, which is slowed to V = 0. 
Q, is such that the loss of momentum is identical in both cases. 

The drag will be 


D = mAV = pQiV 
and the power required to overcome the drag 
HP, = V?/550 


If, now, instead of letting the momentum loss appear in the 
wake we draw off the boundary layer and reaccelerate it to veloc- 


the Editor 


ity V, the power required in the ideal case will be 
HP, = Q, Ap/550 


where Ap = the pressure rise necessary to increase the velocity 
from 0 to V which is precisely (9/2) V?. 


The horsepower required in the second case is, hence, 
HP: = [(0/2)Q,V*] /550 


or one-half the amount needed if the change of momentum is left 
in the air stream. 

This demonstrates the principle involved, and the possibility 
of great power saving, as well as an increase of cruising 
speed. 

Thus, one may state that the theoretical minimum horsepower 
to overcome profile drag is not that which would be needed with 
100 per cent laminar boundary-layer flow but one-half of that 
amount. A suction-type boundary-layer control appears to 
have promising indirect possibilities. 

ALAN Pope 
Associate Professor, Georgia Institute of Technology 
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Chordwise and Beamwise Bending 
Frequencies of Hinged Rotor Blades’ 


GABRIEL HORVAYt 
McDonnell Aircraft Corporation 


SUMMARY 


The &th chordwise natural frequency .\; of a nonrotating blade 
can be written in terms of the fundamental, Nj, as 


Nz = (1) 


The &th frequency WV, of the blade rotating at w r.p.m. can be 
expressed in terms of the nonrotating frequency N; as 


Nz = (Na? + Byw?)'/? (2) 


The significance of these formulas is that A,, B; are rather in- 
sensitive to variations in the mass and moment of inertia distribu- 


_ tions and the plan form of the blade, so that determination of the 


coefficients for one particular blade permits an estimate of the 
various frequencies of other blades of similar design once their 


’ fundamental nonrotating frequency N, is known. 


The coefficients, A;, B, were evaluated for the 23-ft. McDon- 


nell blade. It was found that 
Ao = 0, A, = 1, A, = 3.11 (3a) 
Bo=e, B, =5.06+5.6c, Bs = 13.05 + (3b) 


Here subscript zero refers to the rigid blade mode, and 
e = R,Q/J (4) 


is the “eccentricity parameter’ of the blade. (Ry is the lag 
hinge distance from the center of rotation, Q the static, J the 
inertia moment of the portion of the blade outboard of the lag 
hinge with respect to the lag hinge. For the McDonnell blade 
e = 0.1046.) Calculating the coefficients for beamwise vibra- 
tions and neglecting the flapping hinge eccentricity gives (using 
lower case symbols) 


a=0, a,=1, ae =2.90, a; = 5.90 (5a) 
b=1, 5b, = 5.62, be = 12.6, = 22.4 (5b) 
INTRODUCTION 


Ov ONE OCCASION, when the lag damper of a 23-ft. 
McDonnell blade failed to function properly, 
analysis of the flight record (obtained by a moving- 
picture camera mounted on the swash plate), Fig. 1,1 


Presented at the Rotating Wing Aircraft Session, Sixteenth 
Annual meeting, I.A.S., New York, January 26—29, 1948. 

* This paper is based on an investigation that was carried out 
at the McDonnell Aircraft Corporation under sponsorship of the 
Bureau of Aeronautics, U.S. Navy Department. The writer wishes 
to express his indebtedness and gratitude to C. L. Zakharchenko, 
Chief Engineer of the Helicopter and Jet Propulsion Division of the 
McDonnell Aircraft Corporation, for his never-failing encourage- 
ment, interest, and help, in this and other projects. The writer 
also thanks Jewel Brown Bernard for assistance with the nu- 
merical calculations. 

t Research Engineer. Now Physicist, General Electric Com- 
Pany Engineering General Division. 

t All figures in this paper refer to the 23-ft. McDonnell blade, 
with a lag hinge distance from center of rotation of 16.75 in. 


revealed an appreciable 4w component in the lag vibra- 
tions of the blade. The phenomenon was attributed 
to resonance. Subsequent theoretical analysis con- 
firmed this conjecture. It was found that the first 
chordwise natural frequency of the blade was near 
resonance with the fourth harmonic rotor frequency 
at the particular rotor speed (206 r.p.m.) at which the 
record was taken. 

This result prompted a comprehensive investigation 
of the low-frequency chordwise and beamwise natural 
frequencies of the blade. The values of the frequency 
parameters, A,, B,, a,, b,, quoted in the Summary, are 
the results of this study. 


METHOD OF CALCULATION 


The natural modes of the 23-ft. McDonnell blade 
were determined by the tabular method of Myklestad! 
and Prohl.?** An exposition of this method, as applied 
to rotor blades, was given by the writer in an earlier 
paper.* (The reader is referred to reference 3 for 
further details.) 

Denoting the slope change of the blade at station 
r — / due to a unit moment acting at station r by 


= 1/EI (6a) 


the slope change due to unit shear force, or the deflec- 
tion change due to unit moment by 


Ur—1, Tr = == 1?/2EI (6b) 
and the deflection change due to a unit shear force by 
| = 13/3EI (6c) 


(E is the modulus of elasticity and J is the mean chord- 
wise cross-section moment of inertia between stations 
r — 1 and r), the chordwise vibration amplitudes of 
moment, slope, deflection and shear force at station 
r — l are obtained in terms of those at station r by the 
formulas 


M(r — 1) = M(r) + IS(r) (7a) 
g(r = ¢(r) — rS(r) — 07.1, rM(r) (7b) 
— 1) = yr) — Ie(r — — Sr) 
h,—i,,M(r) 
Sr -)+S(r-) (8) 
where 


** This method can be called the ‘‘moment analysis’ method in 
contrast with other available methods, such as the ‘‘energy”’ 
method of Rayleigh-Prewitt-Wagner-Johnson-Mayne.* 
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Fic. 1. 


Record of test 318-141-1 (hovering), May 26, 1947. Rotor speed = 206 r.p.m., camera speed = 


Recorded points show chordwise deflections of stations 2 and 5 in first cycle. 


184 frames per sec. 


TABLE /, DATA OF THE 23 FT. BLADE. CHORDWISE ViBRATIONS, 


in. 
OF SECTIONS| id ‘ 
IN. iN. FT. | Sevefer.| L8. £8. FT. FT. 108 FT 
276.0\4.0 | O ° ° ° 
276.0 SIF 
269.4 | |.08014 | 72, 136.827) 
262.2 1.725 | 6.1/834/292 | 4.60829I306 | 5021658840 
248.4 1.9 [2.3 |./02/ 2944.80) 3777.14] | 196177 
234.6 2.9 [6.525 70472 | 7.137757 8S4 0.88968 007 
220.5 |“ | 2867.52) 072.763) ./1838/0 
207.0 “  |6.28/0424 | 9,83932050 
193.2 ).7 |" |.09767 | 2917.02) 044.677) .195¥80 
179.4 |6.36097 752 | 7.249965/06 \//.066657899 
165.6 ).6 |“ |.09632 | £77%.09| | ./88637 
151.8 | 7021837026 |/0. 7325/69/38 
438.0 |” |,/000 2984.23 | 1/1F3.6%6 | ./72932 
(24.2 “ S.92947860 | 6. |/0.284726/7 
110.4 |./020 | 2941,9/ |12330.4/0 | ./7512F 
76.6 | 6.$27/3103F | 9.97168629/ 
1.9 |. | 2993.89 |/3228. 559 | ./70//7 
69.0 “ 2¥0 | 6.03//014920 | 9.133 1/9614 
FIR | |/4270.919 | 154423 
41.4 S96 592335F | 8.990399 
27.6 A304 | 08S. 99 |/$379.518 | 
16.75 |2,00304Y/47 | | 5456377088 
/6.75|.0607, |. 8060 ° | 153872 
R 
‘ elop were it to act at the blade tip R, and 
R 
is the inertia shear, and C(r) = Lirr/R (10c) 
r 
R 
Se(r) = —C(r)e(r) + > Ty/R (9b) is the centrifugal tension at station r.* 
is the centrifugal shear. Here The determination of S, Eq. (8), is carried out more 
—* (10a) conveniently in the form 


is the ratio of vibration frequency to rotor frequency, 
t(r) = 


is the centrifugal force of the blade segment d (of mass 
density, m’, lumped into station r) which it would de- 


(10b) 


S= 


where 


Si +. S2 


(11) 


* Evidently, the use of the above formulas implies the assump- 
tion that all deflections are small so that their higher powers are 


negligible. 
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BENDING FREQUENCIES OF HINGED ROTOR BLADES 


TABLE 2. DETERMINATION OF CHORDWISE 


BLADE MODE. =223RPM 


Paton ~ 
B70” Coefficien7 of Coefficient of fomer- 
Ltem Y= sea? 
4 Oo 23.0 40 
x2.975° 
Mm 
40 
-.575 23.0 
—288-9730 U557.20 
S, | — M26. 827 ° 
S | 4/5. 757 NSS7.20 
M | ~2442.1/8082S5'| 17736-1709 
¥ 4006524224 | --05325896740 
4 ~2-304/44838 |23.03383540 +652 
Sz |-3237.036834 | 4/048. 62108 
S, |— 380/.789953 201.16 
S | 7040. 826787 | 4/249. 68803 
| 48/2004.0/3 
Y | | ~S0./22/87/7 
| 2/97.3802 345° —-.333 
S, |-283249.482j | 764202-357/ 
S, |-/97724.786/ | 7/5287. 6732 
S | 682 | 
! 
M |-243166.87/ | 8,674835. 682 
Y | 23-4406738E | —84.54097653 
4 |-/06.8289388 | 370. 4783857 
Sz |- 778/64. 5026 | 2,750002. 9765" 
S, |-360506.2659 | 1300/77. 778 
S |-¥58670.769 | 4, 050202.743 
=.0607 
| 12,3362/8-F6 |8/57-8 
Y | 29.3S5944 -105.5 8403/5 —-/326 
~130- 5365365 | 455. 8824103 oO 


Si(r) = —C(r)e(r) 


R 
S:(r) = (1 + n?) ry/R 


(12a) 


Letting J(r) represent the beamwise cross-section 
moment of inertia of the blade and writing 


= —C(r)¢(r) 
R (12b) 
S:(r) = 
the formulas (7), (11), and (12b) are those for beamwise 
vibrations. 
At the blade tip, the moment and the shear vanish,* 
Mr = Sr = 0 (13) 
The tip deflection 
can be assumed to be unity, while the tip slope 
(14b) 
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must remain an unspecified quantity. Carrying out 
the calculations, station by station, from tip to root, by 
formulas (7), (11), and (12), for an assumed vibration 
frequency vy = nw, the condition of no hinge deflection 


Ja (15) 


fixes the value of the tip slope ® and also of the moment 
— My, which must be applied at the hinge to provide 
blade oscillations of unit amplitude. Resonance ob- 
tains at the frequencies v for which 1/, = 0. 

Thus, in the example of Table 2 (this is based on the 
chordwise blade constants of Table 1), it is found that 
for the vibration frequency »y = 3w = 669 cycles per 
min., tip slope ®, and tip deflection RY = 1 ft., the 
condition of vanishing hinge deflection 


Ya = —130.53653656 + (455.8824103/23) = 0 
(16) 

fixes the amplitude of the tip slope at 
gr = ® = 0.1518423492 (17a) 


and hence the hinge slope at 


= —0.1326 (17b) 
and the hinge moment at 
My = 8,157.8 Ib.ft. (17c) 


A calculation of My vs. v indicates (Fig. 2) that 
resonance obtains at the frequencies v (cycles per min.) 
equals 


M = 0.3234 X 223, N, = 


~10, 


7 
| 
| 
= | \ | 
/ 
T 
\\ / / 
\ / 
| | \ \ / 
© CALCULATED POINT 7 
——— BLADE STATIONARY 
—— BLADE ROTATING AT 223 RPM / 
-30, 
\ / 


Fic. 2. Chordwise frequency response curves. My cos 
vt is the moment that must be applied at the hinge to maintain 


* Subscript T refers to blade tip; subscript H, to hinge values. tip vibrations of 1-ft. amplitude. 


= 
| 
3.64 223, 
Nz = 9.38 X 223 (18a) 
| 
| | | 
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Fic. 3. Beamwise frequency response curves. 


for a blade rotating at w = 223 r.p.m., and at the fre- 
quencies 


No = 0, M = 2.76 X 223, Nz = 8.57 X 223 (18b) 


for a nonrotating blade. 


As an example to beamwise frequency calculations 
the reader is referred to the columns ; & 4, and Y; & 
Y, of Table 4a, reference 3, which indicate that, for a 
tip slope ®, tip deflection RY = 1 ft., and vibration 
frequency vy = 2w = 446 cycles per min., the condition 
of vanishing hinge deflection 


Ya = —3,159.091619% + (9,482.335974/23) = 0 (19) 


implies a tip slope 


® = 0.1305044382 (20a) 
and a hinge moment of 
My = 4,378.3 Ib.ft. (20b) 
In this case one finds (see Fig. 3) that 
io = 223, v1 = 2.54 X (21a) 
ve = 4.42 X 223, v3 = 7.16 X 223 


are the rotating, and 


= 0, = 0.91 X 223 
= 2.64 X 223, = 5.37 X 
are the nonrotating natural frequencies (cycles per 


minute) of the blade. 


1948 


CHORDWISE BLADE MopEs 


The blade deflection curves, y(r), at the first two 
natural frequencies are shown in Fig. 4. It is seen that 
the mode shapes of the rotating and nonrotating blade 
differ but slightly, and the nodes are practically coin. 
cident. 

The calculated first natural frequency 


Ni cate. = 2.76 X 223 = 615 cycles per min. (22a) 


T 
‘——— BLADE STATIONARY | 
—— BLADE ROTATING AT 223 RPM 


| 


DEFLECTION 
/ 
\ 


(a) 
| | 
‘ 
| 
z J —— BLADE ROTATING AT 223 
61 A \ 
2 2 4 6 8 r 
(b) 


Fic. 4, Chordwise natural modes. (a) Fundamental mode, (b) 
second mode. 
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Fic. 5. Variation of the first chordwise natural frequency with 
- the speed of rotation. 
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differs from the experimental value 


Niexp. = 650 cycles per min. (22b) 


indicating that the chordwise stiffness of the blade (EJ 
of Table 1) has been underestimated in the calculations. 
Since the present paper aims primarily at determina- 
tion of frequency ratios rather than of absolute fre- 
quencies, this discrepancy is of minor import. By 
Eqs. (1) and (18b) one obtains the frequency ratios 


Ao = 0, A; = i. A, = (23) 


for the nonrotating blade, and by Eqs. (2) and (18a) the 
correction factors 


= V1 + (Bywo?/N 2) 
By = 0.1046, B, = 5.63, By, = 14.45 


for the rotating blade. 

To test the validity of the form of formula (24a), 
My was calculated at w = 111.5 r.p.m. for an assumed 
vibration frequency 


= (615? + (5.63 111.5%)]'”* = 669 cycles per min. 


The small value, Wy = 30.4 lb.ft., obtained for the re- 
sidual moment, confirms the correctness of formula 
(24a).4 

Fig. 5 illustrates the dependence of M. on w. It 
shows that, for a nonrotating frequency of Ni = 650 
cycles per min., resonance with 4w occurs at 203 r.p.m. 
(This compares with the 206 r.p.m. at which the flight 
record, Fig. 1, was obtained.) 

In order to determine the dependence of the coef- 
ficients B, on the blade eccentricity ¢«, the frequency 
response curve of the rotating blade was recalculated 
for lag hinge separations R, = 0, 33.5, and 67 in. The 
calculations led to the consistent set of values [Eq. 
(3b)] quoted in the Summary.* The value 


De = 


(24a) 
(24b) 


(25) 


is the well-known Coleman formula, which plays an im- 
portant role in ground vibrations. 

It is well to point out at this place that the analysis 
presented above implies that no coupling exists be- 
tween lag motion and blade pitch angle. It was shown 
by the writer’s colleague, A. Gail (unpublished result), 
that in case of coupling the simple Coleman formula 
must be replaced by 


(26) 


where Q, is the aerodynamic torque; #, the blade 
pitch; and ¢, the lag angle.** Correspondingly, it 


* Evidently, two RX, values are sufficient to determine the two 
terms of a linearly assumed B,. The extra two values served to 


verify the correctness of this assumption. 

** An expansion of Q,, Eq. (75) of reference 3, in powers of ¢, 
and transfer of the linear term to the left-hand side yield immedi- 
ately Gail’s formula, Eq. (26). 
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Fic. 6. Angle ¢# of chordwise hinge deflection. 


BENDING MOMENT AT HINGE, LB FT. 


BENDING MOMENT AT HINGE, LB. FT 


10,000} 


10x 223 


Chordwise frequency response 1 of damped blade; 

D = 400 lb.ft. per rad. per sec. Also undamped moment My, 

and damper moment Mz, curves. (The signs of My and Mg are 

disregarded.) (a) Blade stationary, (b) blade rotating at 223 

r.p.m. 


Fic. 7. 


must be expected that the higher B, will also be in- 
fluenced by the slope of the #(¢) relation. [For the 
McDonnell blade #(¢) is essentially a constant. } 


ACTION OF THE BLADE DAMPER 


In what has preceded it was tacitly assumed that no 
damper acted at the lag hinge. If there is a viscous 
damper attachment of damping constant D, then 


(27) 


represents the moment absorbed by the damper. The 
factor j indicates that D¢y is 90° out of phase with the 
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' = -— The hinge slope ¢y is plotted vs. v in Fig. 6; © and, 
for comparison, also and M,, with signs disregarded Th 
> a | are plotted in Fig. 7a (nonrotation) and Fig. 7b (rota. ip 
8 P = a tion at 223 r.p.m.) for a damping constant of D = 4 pot 
Ib.ft. per rad. per sec. tion. 
OR MN is always equal to, or larger than, My and \, =. 


Fic. 8. Rigid rotating blade at Points of tangency occur only where gy or My vanish, | ihe 
ee P It is noted for the rotating blade that up to about p = , | bendi 


the damper absorbs little vibration. N evertheless, detent 
| | | | this small damping plays an important role in prevent. }, = on 
ee on oe ing ground vibrations. At the higher resonance fre. } load, 
5 | quencies is extremely effective in limiting the vibra. youd 
Bo | {| | uR tion amplitudes to small values. It is noted that the ie 
1 resonance frequencies (9% = minimum) are almogt 10 chi 
i= ee entirely unaffected by the presence of the damper rey 
| | (a) | (For exceedingly large damping constants, not achiey. theory 
able in practice, this would no longer be true. tion 0 
bl pract th Id no longer be true.) The 
same is found true also for the deflection curves. - te 
| | Prof 
| | | | | BEAMWISE NATURAL MODES modul 
| In Fig. 8 the rigid mode, = and the static de- | a’bett 
flection curve, = 0, are shown for the rotating blade, | modul 
7 7 The first, second, and thiid beamwise natural modes, 
| | nonrotating and rotating at 223 r.p.m., are shown in | furthes 
| | Figs. 9a-9c. The discrepancy between the rotating | ‘sinvo 
| | | and the nonrotating modes seems to increase with the wall” 
| A somewhat different conclusion is reached by Simpkinson, Fron 
Eatherton and Millenson,' for cantilever blades. paper, 
wr, difficul 
é REFERENCES fusion 
zybski’ 
1 Myklestad, N. O., Vibration Analysis; McGraw Hill Book } an actt 
= 1 Company, Inc., New York, 1944. tion th 
— ey | 2 Prohl, M. A., A General Method for Calculating Critical Speed | centrall 
= Bon ng of Flexible Rotors, Journal of Applied Mechanics, Vol. 12, p. A-142, | °rder 0 
4 3 Horvay, G., Stress Analysis of Rotor Blades, Journal of the — 
column 


Fic. 9. Beamwise natural modes. (a) Fundamental mode, (b) Aeronautical Sciences, Vol. 14, p. 315, 1947. Errata to the above | ip.+ in 


second mode, (c) third mode. paper, ibid., Vol. 14, p. 544, 1947.7 it is of 

4 Den Hartog, J. P., Mechanical Vibrations, p. 310; McGraw | we coul 

undamped hinge moment My. The moment that Hill Book Company, Inc., New York, 1944. of a fin 
t be supplied at the hinge to shake the blade ti 5 Simpkinson, Eatherton, Millenson, Effect of Centrifugal | an enti 
ore a PP : - 8 : Shake the Diade Up Force on the Elastic Curve of a Vibrating Cantilever Beam, N.A.C.A. | attacke 
with unit amplitude is —Wy, — jM,,and hasthe mag- TN, 1204, 1947. his 190: 
nitude : small ec 
+ In correction (15) of this reference, (§)4, which is the coefli- | because 

orders o 


= Vu a+ (28) cient of the sin 2y term of §, should read —2wAo. 


Letter to the Editor 


Dear Sir: loads would be much too high. It is unfortunate that the coefl- 
cient of fixity was ever introduced into the column literature. 
The following notes were suggested by Professor Wang’s paper _ Since end restraint reduces the effective column length, its effect 
“Inelastic Column Theories and an Analysis of Experimental should be taken care of by a factor that operates on the colum 
Observations” (JOURNAL OF THE AERONAUTICAL SCIENCES, May, length, not on the load or stress. Such a factor has been used by 
: 1948). many authors, in the form, k = Lo/L. I am heartily in favor of 
Fs In Eqs. (16), (18), and (19) the ‘‘coefficient of fixity” was used dropping C entirely and substituting an effective length factor. 
‘ as a multiplying factor for the critical column load in the inelastic | (Incidentally, I can claim priority over Wang on this error, hav- 
range. If C were actually used in this manner, the predicted ing done the same thing in my book, Basic Structures!) 
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The second point concerns the use of the term ‘“‘neutral axis.” 
In his 1909 thesis von Karman pointed out (in a footnote) that 
the neutral axis as he used it was not “‘stress free,’’ since axial 
compressive stresses were also acting over the entire cross sec- 
tion. If the curve of actual stress distribution is known, the 
“geutral axis’ is found at the point of intersection of this curve 
with a horizontal line representing the average axial stress. For 
the elastic case this coincides with the neutral axis for pure 
bending, but in the inelastic case it does not; neither does it 
determine the location of the resultant axial force for a constant 
stress distribution (this remains at the centroid). If the bending 
is permitted to proceed simultaneously with the increase of axial 
load, the “neutral axis’’ will not be stationary but will move 
away from the centroid as the loading increases. The only 
physical meaning that could be attached would be that it de- 
termines a line along the length of the column for which there is 
no change of axial stress. But this requires the assumption of a 
constant average axial stress during bending, the very assumption 
that is now agreed to be wrong and on which the double-modulus 
theory was based. It would seem desirable to restrict the defini- 
tion of the term ‘‘neutral axis’? to pure bending. This would 
force realization that its use in combined bending and axial load- 
ing involves treating the bending separately. 

Professor Wang has attempted to improve the simple tangent 
modulus formula by introducing two modulii, one higher and 
one lower than the tangent modulus. This is done to obtain 
a’better approximation to the true stress distribution. These 
modulii are then used in the original double-modulus formula, 
again introducing the original assumption of a neutral axis. 
Furthermore, the use of two modulii aE; and bE; requires a 
further assumption as to how much differential strain (bending) 
isinvolved. If we use the idealized column as a basis for analysis, 
we must assume, as Wang has assumed, an “‘infinitesimally 
small” disturbance. This implies that the differential strain 
approaches zero as a limit, bringing us back to the tangent 
modulus again. 

From’ the correspondence that followed the publication of my 
paper, “Inelastic Column Theory,’ I can see that much of the 
difficulty associated with buckling theories comes from a con- 
fusion of ‘‘orders of abstraction” (to borrow a term from Kor- 
zybski’s Science and Sanity). Since it is impossible to analyze 
an actual column, we must create a mental picture or abstrac- 
tion that may be more or less idealized. The so-called ideal, 
centrally loaded column in the elastic range represents a high 
order of abstraction and leads to the original Euler formula. A 
closer approach to reality is-obtained by replacing the elastic 
modulus by the tangent modulus, but the idea of the ‘“‘perfect”’ 
column is still retained. The degree of eccentricity or disturbance 
that must be assumed in such an analysis is truly infinitesimal; 
it is of the order of the disturbance needed to upset a needle if 
we could succeed in making it stand on end! The introduction 
of a finite eccentricity or ‘‘disturbance’’ places the problem on 
an entirely different level of abstraction. It must then be 
attacked by methods similar to those used by von Karman in’ 
his 1909 paper. Most of the arguments about the effects of 
small eccentricities on the ‘“‘column theory”’ can never be resolved 
because they are not arguments at all, simply confusions of the 
orders of abstraction involved. 
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Reduced Interaction Curves for Combined Bending and Axial Loading 


Fig. 15 


I do not believe that any additional theories or modifications 
are needed to explain the fact that test values often exceed the 
values predicted by the tangent modulus theory in the region 
above the knee of the stress-strain diagram. The remarks 
about Fig. 10, in my original paper, together with Dr. von 
Karman’s discussion, should mdke this quite clear. By the 
same reasoning, it can be seen that, in the region of the knee of 
the stress-strain diagram, where the tangent modulus decreases 
rapidly with increasing strain, the inevitable initial eccentricities 
and imperfections will have their largest effect and test data 
will tend to fall below the tangent modulus values. 

Finally, I should like to suggest that the problem of finite 
eccentricities be attacked by an entirely different engineering 
procedure, the use of the interaction curve for combined axial 
loading and bending. In a recent A.S.C.E. paper (‘‘Applied 
Column Theory,” not yet published) I showed how this could be 
done by adding lines of constant eccentricity ratio to an inter- 
action chart. This method covers the entire range from pure 
axial loading to pure bending, is dimensionless, and can be re- 
duced to a simple form as shown herein by Fig. 15 from the 
above paper. This chart is based on N.A.C.A. Technical Note 
No. 307, in which the tests were made by transverse loading at 
the third points. For eccentric columns (bending moment ap- 
plied at ends) ‘the curves would show more sag. I am now at- 
tempting to collect test data on eccentric columns for which the 
bending modulus of rupture has also been determined and would 
appreciate receiving any information of this type. 

F. R. SHANLEY 
Consultant on Structures 
Los Angeles, Calif. 
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